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Chapter 1 

Introduction to Hydrodynamics 


Sangyong Jeon^ and Ulrich Heinz^ 

^Department of Physics, McGill University, Montreal, QC, Canada 
^Department of Physics, The Ohio State University, Columbus, Ohio, USA 


1. Introduction 

Application of hydrodynamics in high-energy physics has a long and illustrious his¬ 
tory starting from L.D. Landau’s seminal paperl^^ In its history of more than half a 
century, many papers have been written on a broad spectrum of topics, too numer¬ 
ous to list them all here. In this review, our emphasis will be more on the basics 
of the theory of hydrodynamics than to report on the current phenomenological 
status, of which several excellent reviews already exist (for instance, see Refs.l^. 

Recent ultra-relativistic heavy ion collision experiments at the Super Proton 
Synchrotron (SPS), Relativistic Heavy Ion Collider (RHIC) and the Large Hadron 
Collider (LHC) have demonstrated beyond any doubt that Quark-Gluon Plasma 
(QGP) is being created in these collisions. Unfortunately, direct access to the QGP 
properties such as the temperature, equation of state, transport coefficients, etc. is 
not very feasible. The only experimentally accessible information is contained in the 
spectra of the final state particles. To connect them to the QGP properties such as 
above, one must use theoretical models. It would be wonderful to have an analytic 
or a numerical method that can calculate evolution of heavy ion collisions from 
first principles. But this microscopic, non-equilibrium, many-body QGD problem 
is currently intractable. What is tractable is the coarse-grained collective motion 
of the system as a fluid after the local thermal equilibrium is established. Since 
the properties of QGP we are after are mostly (local) equilibrium properties, it is 
natural that the dynamics of collective motion - hydrodynamics - is an integral 
part of the theoretical modelling. 

What has been exciting and interesting in QGP research is the close discourse be¬ 
tween the experiment and theory. In elementary particle experiments, perturbative 
QGD is being tested with amazing successes. More and more precise perturbative 
QGD calculations prove to describe experimental data more and more accurately. 
In contrast, QGP research is much more dynamic. For instance, before the dis¬ 
covery of QGP, theoretical expectation was that QGP would be a weakly coupled 
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plasma of quarks and gluons based partially on the fact that the QGP properties 
seem to approach about 80 % of the Stefan-Boltzmann limit rather quickly on lat¬ 
tice, around 2T(W But almost from the day-1 of RHIC operation, strong elliptic 
flow quickly proved this initial expectation not very viable. QGP around the tran¬ 
sition temperature turned out to be the most strongly coupled many-body system 
ever observed. Soon after, the authors of Ref.l^ used string theory techniques to 
calculate the infinite coupling limit of the shear viscosity and came up with a sur¬ 
prising result that the limit is small, but has a non-zero lower bound. Subsequent 
Hydrodynamic calculations then demonstrated the importance of small but finite 
shear viscosity in understanding the RHIG flow data.^^^ Gomparing the ensuing 
LHG predictions with the LHC data now confirmed the expectation that as the 
temperature increases, shear viscosity of QGP should also increase.^^^^^ 

All these connections between exciting theoretical developments and experi¬ 
ments cannot be made without hydrodynamics. More recently, the systems created 
in the highest multiplicity proton-proton collisions and proton-nucleus collisions 
were also seen to exhibit strong collective behavior.^i^®^ This is deeply puzzling as 
the size of the system ought to be too small to behave collectively. It is hoped that 
more thorough investigation of the possible origin of the collectivity in such small 
systems can illuminate the inner workings of the QGP formation greatly.^^^ 

As mentioned in the beginning, the aim of this review is the introduction of 
the theory of the hydrodynamics in ultra-relativistic heavy ion collisions. This 
actually entails a large number of disciplines in addition to the relativistic fluid 
dynamics. Our plan for this review is as follows. Section 2 contains the basic 
concepts of hydrodynamics and their definitions. In sections 3, second order viscous 
hydrodynamics is derived from a very general linearly response theory of conserved 
currents. Section 4 discusses how coarse-graining of kinetic theory can result in more 
general form of viscous hydrodynamics. In section 5, various numerical techniques 
needed to implement relativistic viscous hydrodynamics in ultra-relativistic heavy 
ion collisions are discussed. We conclude in section 6. 

2. Hydrodynamic Form of the Stress Energy Tensor and the Net 
Baryon Number Current 

Hydrodynamics is all about flow of conserved quantities. In this review, we strictly 
deal with relativistic hydrodynamics. Therefore, unlike the non-relativistic case, 
mass is a part of the energy budget. In the Minkowski coordinates, conservation 
laws in their local form are 

= 0 

= 0 ( 1 ) 

where is the stress-energy tensor and the roman letter i on the current labels 
any other conserved charges such as the net baryon number, net electric charge, etc. 
For the bulk evolution in relativistic heavy ion collisions, usually only the net baryon 
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number current, is considered. If needed, additional electric current and the 
strangeness current can be easily accommodated. 

Using the divergence theorem, the integral form of the conservation laws read 



( 2 ) 


where in the right hand side the integration is over the boundary of the volume V 
assuming that the size and the shape of the volume is independent of time. 

This form admits a very physical interpretation that the rate of change of the 
conserved quantity in a fixed volume equals the net current entering the volume. 
Hence, the dynamics of conserved quantities are governed by the dynamics that 
governs the currents. In essence, hydrodynamics is all about the dynamics of the 
currents. 

Hydrodynamics is useful because it is a coarse-grained theory. When a system 
contains too many particles, it becomes difficult to follow microscopic details of the 
system. When the system contains sufficiently many particles, the system again 
starts to admit analytic studies because thermodynamic concepts start to apply, in 
particular the static equilibrium. A system in static equilibrium is characterized 
by only few quantities such as the temperature, collective velocity and chemical 
potential. These control the energy density, momentum density and charge density, 
respectively. The price to pay for this simplification is that questions on short time 
scale phenomena or short length scale phenomena cannot be answered any longer. 

The systems we would like to study, the ultra-relativistic heavy ion collisions, 
do contain a large number of particles, but they certainly are not static. In fact, 
they will never actually reach the state of static equilibrium. Nevertheless, if one 
is interested only in the coarse-grained collective motion of the system, the concept 
of local equilibrium may still apply provided that the expansion rate is much slower 
than the microscopic interaction rate. If one considers a macroscopically small but 
microscopically large fluid cell around a position x at a given time t, then within a 
macroscopically short but microscopically long time scale, the time averages should 
approach the static equilibrium values according to the ergodicity hypothesis of 
statistical mechanics. More details on the length and the time scale analysis can be 
found in Section lUTI 

When the local equilibrium is reached, it becomes meaningful to describe the 
system with the local temperature T{t, x), the collective velocity u^{t, x) of the fluid 
cell and the local chemical potential fasit, x). One can then study dynamics of only 
those few thermodynamic quantities. Since are basically the Lagrange 

multipliers to fix the average energy, momentum and net charge, it is natural that 
we turn to the conservation laws for their dynamics. The goal of hydrodynamics is 
then to study collective motion of a system using the conservation laws with only 
the statistical inputs from the underlying microscopic theory. 
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In an dynamically evolving fluid, the concept of the local rest frame is essential 
in order to apply the concept of local equilibrium. However defining the collective 
velocity of a fluid cell turned out to be quite intricate. This is not so simple even 
for the simplest system composed of a single kind of non-interacting particles in 
the non-relativistic setting. One kind of collective velocity comes from the average 
momentum 


Up = 


EiW 


( 3 ) 


or from the mass current. Here the sum is over all particles in the given fluid cell. 
Another comes from the energy-weighted average 


Ue 


Ez(wvf/2) 


( 4 ) 


or from the energy current. If the particles in the system have an additional con¬ 
served charge, B, and the net charge is non-zero, then one can define yet another 
collective velocity by performing the charge-weighted average 


ub 




( 5 ) 


Here Bi is the conserved charge of the i-th particle. 

For a state in static equilibrium, all three collective velocities above coincide 
because they must all vanish. However, they do not necessarily coincide to an 
observer moving with a uniform speed — uq. The mass weighted average velocity 
and the charge weighted average velocity are both uq- But the energy weighted 
average velocity 


, X^^(m(v, + uof/2)(v, + uo) 

' E,m(v, + uo)^/2' 


^ Uo 


( 6 ) 


clearly does not coincide with uq. There is no unambiguous choice of the flow 
velocity even in this simple case. One must choose among these options what will 
be regarded as the flow velocity. 

In the relativistic setting, mass is a part of the energy. Hence, there are only 
two options for choosing the flow velocity: The energy current or the net baryon 
current. One must choose one of these velocity options in order to decompose T^'' 
and into a useful hydrodynamic form. The net baryon number is relatively small 
in ultra-relativistic heavy ion collisions. Furthermore, the flow observables we are 
interested in are mostly patterns in energy and momentum distributions. Therefore, 
there is no real benefit to choose the net baryon number collective velocity as long 
as the heavy ion analysis is concerned. 

Choosing to follow the energy current the flow velocity for the energy current 
is defined by the eigenvalue problem 


= evF- 


( 7 ) 


^This choice of frame is often referred to as the Landau-Lifshitz frame. If one choose to follow the 
charge current, it is referred to as the Eckart frame. 
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where = T^°‘gai, and the flow vector is normalized to = 1. 

We use 5 ^ 1 / = diag(l, —1, —1, —1) throughout this review. Note that while is 
a symmetric matrix, is no longer a symmetric matrix. Therefore, there is no 
guarantee that the eigenvalues are real. But any physically consistent system should 
admit a positive eigenvalue £ and the associated time-like real eigenvector . 
Decomposing using e and one gets 

= eu^^u'^ + \ (T“ - £) (8) 

where the local 3-metric is given by 

= gt^'' - (9) 

The residual shear tensor is symmetric transverse = 0 and 

traceless = 0. Hence altogether the expression ^ has the required 10 

independent degrees of freedom; the local energy density e, the local fluid velocity 
u, the trace T°^ and the residual shear tensor The net-baryon current is 

= + ( 10 ) 

where Jg = pbu^ is the ideal fluid part of the current and Vg is a space-like vector 
satisfying the transversality condition = 0. It has the required 4 independent 

degrees of freedom; the local net baryon density ps and the residual vector V^. 

So far we did not use any thermodynamic information. We will do so now 
to re-write the trace T°^ in a more physical form. A static medium at rest has 
'^eq = diag(e, P, P, P) where P is the pressure. Therefore, the trace should contain 
the equilibrium piece Pf^uT^q = s — 3P. Furthermore, the thermodynamic identities 

dP = sdT + psdpB (11) 

ds = ^de-^dpB (12) 

where s is the entropy density, indicate that the pressure P is a function of the 
temperature T and the baryon chemical potential pB, and they are in turn function 
of the energy density and the net baryon density. Hence we must be able to find P 
as a function of £ and pB- 

P = P{e,PB) (13) 

This relationship is known as the equation of state. Writing = £ — 3(P -I- H), 
the stress-energy tensor then becomes 

_ nA'"'^ -h TT^^’' (14) 

where = eu^u'^ — P{e, pb)^^'' is the ideal fluid stress-energy tensor. 

From the arguments presented above, it should be clear that the residual scalar 
term H and the tensor term as well as the vector in the baryon current must 
vanish in the static equilibrium limit. As these quantities represent deviation from 
equilibrium, the size of these terms will depend on how fast the local equilibrium is 
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achieved. If local equilibration is instantaneous in the macroscopic time scale, then 
fluid cells will always be in strict local equilibrium and hence 11 = = 0. 

Microscopically, this would happen if scattering cross-section is large so that the 
mean free path is much shorter than any macroscopic time scale or length scale 


(c.f. section 4.1). If this is the case, the number of unknowns (e, u,/9b) matches 
the 5 conservation laws and one has the ideal hydrodynamics whose dynamics is 
completely specified by 


d^ieuV -P{e,pB)A>^n = 0 


(15) 


and 


df,[pBU^) = 0 (16) 

The information on the underlying system is only in the equation of state P{e, pb)- 


3. Hydrodynamics from Linear Response Theory 
3.1. Linear Response Theory 

In realistic systems, the approach to the local equilibrium is never going to be 
infinitely fast. Therefore, 11, and Vg cannot simply be set to vanish, although 
if the system is conformal one strictly has 11 = 0. This is because the trace must 
vanish = 0 in a conformal system. When any of these quantities are non¬ 
zero, the system is out of equilibrium. Therefore, local entropy must increase. The 
evolution equations of these quantities are then necessarily of the dissipative type. 
In this and following sections, we use linear response theory to obtain such equations 
for dissipative hydrodynamics. In section]^ kinetic theory approaches that can go 
beyond the near-equilibrium restriction of the linear response theory is discussed. 

To gain more insights on the behavior of the dissipative quantities, we start 
with a system very slightly out of equilibrium at t = 0. We then consider how the 
system approaches the equilibrium. This is the realm of the linear response theory. 
Full analysis of the quantum linear response theory can be found in any number of 
standard text books (for example, see Ref.^^. Here we will only go over the main 
ideas. 

Suppose that at the remote past, t = —T, the density operator had the equilib¬ 
rium form pq = e~^^°jZo where Hq is the system Hamiltonian and Zq = Tre“^^° 
is the partition function. Then a force term is adiabatically turned on f(t,x) = 
at an infinitesimally slow rate e. At t = 0, the force term is turned 
off and at this point the system is out of equilibrium. The full time-dependent 
Hamiltonian for this process is 

H{t) = Hq - J (fxA{x)f{t,x.) (17) 

Here A represents a set of Hermitian operators for the conserved quantities we are 
interested in. Namely, and pB- Hence, the expressions presented below are in 
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general matrix expressions. 

Treating 5H = —f d^x f(t,x)A(x) as a perturbation, the formal first order 
solution of the quantum Liouville equation idtp(t) = [H{t),p(t)] is given by 

Spnit) =-i J dt'[SHH{t'),po{0)] (18) 

where Sp{t) = p{t) — po and the subscript H denotes Heisenberg picture operators 


dH{t) = 


(19) 


Using 5H from Eq.( 17), the deviation of an observable A from the equilibrium value 
is found to be 


5{A[t,x)) = J d'^x'G^^{t — t',x —x!)0{—t')e'^*' f{x') (20) 

where t > 0 and 

(i - , X - x') = (t - t') Tr (|po [Iff (t, x) , H// (t', x' )] ^ (21) 

is the retarded response function. We also took the —T —>• —oo limit. 

Suppose that one can find an operator Da for which is the generalized 

retarded Green function, 


DAGR^{t — t',x — x') = dAS{t — t')S{x — x') 


( 22 ) 


where dA can cont ain a finite number of derivatives. For t > Q, t and t' can never 
be the same in Eq.(20). Hence 5 {A) satisfies the evolution equation 


DAS{A{t,x)) = 0 


(23) 


for t > 0. Therefore finding the pole structure of the response function is equiva¬ 
lent to finding the evolution equationEnMl We will use this to find hydrodynamic 
equations in the following sections. 

For further analysis, it is useful to define the spectral density 

p^^iu;, k) = I d^x ([iff (t, x), it^(O)]) (24) 

Using the thermal average ((• • •)) = Tr/5o(- • •), it is not hard to show (for instance, 
see Ref.li^ 

p^^(w,k) = {2TT)‘^S{k - Pm + Pn) 

Zjq 

m,n 

where k = (oj, k) and jm) is the simultaneous eigenstate of the system Hamiltonian 
Hq and the total momentum P of the system with the eigenvalue Pm = {Em-, Pm)- 
From this expression, we can derive 

p^^{-uj, -k) = -p-^^{uj, k) (26) 



by exchanging m and n. When the underlying equilibrium system is isotropic, then 
p'^'^(a;,k) must be a function of jkj only. Hence the spectral density is an odd 
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function of w, w, k) = We will use this property often in later 

sections to parametrize the analytic structure of the response functions. 

In terms of the retarded correlator is 



^ p^^(cc',k) 

27r Lo' — uj — ie 


(27) 


The imaginary part of the retarded correlator directly gives the spectral density 




(28) 


It is also useful to know that the Euclidean correlator is given by (for derivations, 
see, for instance, Ref.Eal) 




2tT Uj' — iuJn 


(29) 


where = 2'KnT is the Matsubara frequency. One important fact we will often 
use in the following sections is that the a; —>■ 0 limit 


G^^(0,k) = G^^(0,k) 

r doj'p^^{uj,k) 

I 271 Uj' 




g-HB, _ 5-SE. 


1 - 2 
(27r)^<5(k - Pm + p„) (n|A|TO) 




(30) 


is real and positive and function only of the magnitude of k0 


3.2. Baryon Density Diffusion 

From the previous section, it is clear that the analytic structure of the response 
function determines the evolution of small disturbances. In this and the following 
sections, we show that this fact combined with the conservation laws is powerful 
enough to produce dissipative hydrodynamic equations. 

The analysis in this section and those up to section [STS] closely follow the unpub¬ 
lished note by Laurence G. Yaffe (private communication, see also RefsBHlHl) in a 
simplified form. Additional discussion on the 2nd-order formulation of dissipative 
hydrodynamics is given in section |3.6[ 

We start with the net baryon number conservation which is the simplest to 
examine. Suppose we set up a system where only the net baryon density (pb(0,x)) 
is non-uniform at f = 0. The perturbing Hamiltonian is 

SH{t) = - J d^x pb (x) e"* pb (x) 

^To see this, just exchange m and n. 


(31) 
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which results in the following linear response in the mixed space of t and the 
wave vector k, 


/ OO 

dt' G™ {t-t', k) 

-OO 


(32) 


for t > 0. 

Applying the current conservation, dtps = ~V-Jb, to the retarded correlation 
functions G^'^(t,x) = x), J^(0)]) results in the following relationships 

between them in the frequency-wavevector space 


wG?j°(a;,k) = fc,G?j*(cu,k) 
ujG°^{uj,k) = fciG)((w,k) 


(33) 

(34) 


provided that [J)^(0, x), J^(0, x')] = 0. The underlying isotropic equilibrium per¬ 
mits decomposition into transverse and the longitudinal parts 


G%{uj,k) = k^yGL{oj,k^)+S^^GT{uj,k) 


(35) 


where k = k/|k| is the unit vector and jb = Jb _ jg ^110 transverse projector. 

(36) 


Combining Eqs.(33) and (34) gives 

uj^G°]^{Lu,k) = k‘^GLiuj,k) 


What we are interested in is the behaviour of GL(a;,k). 


Consider the small uj limit of Gl first. From Eq.(30), we know that 
5oo(k)=G™(0,k) = G™(0,k) 

is real and positive. The small oj limit of Gl can be then expressed as 

w^5oo(k) 


GL(cc,k) 


k2 


(37) 


(38) 


Now consider taking the k —0 limit with a fixed ui ^ 0. The retarded density- 
density correlation function in this limit is 


G?j°(a;,0)=i J dte J d^x {[pB{t,x), pb{0)]) 


= i 


dte ([Qb,Pb(0)]) = 0 


(39) 


where we used the fact that Qb = f d^x pB(t,x:) is the net baryon number operator. 
Since the net baryon number is conserved, Qb is independent of time. In particular, 
it can be evaluated at t = 0. Therefore, the commutator vanishes. This then indi¬ 
cates that GB(a;,k) is well behaved in the zero |k| limit with oj ^ 0. Consequently, 
the OJ 0 limit and the k —> 0 limit do not commute, indicating the presence of a 
massless pole. We also know that the imaginary part of G^j* (the spectral density 
p^^) has to be an odd function of oj since isotropy in space demands that it be a 


function only of k^ (c.f. Eq.(26l). 
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The most general form of Gl consistent with the above conditions is 
. a;2(goo(k)+iwA(a;,k)) 

K) — 


(40) 


— iuj/D[iu^ k) — k) 

The functions D{u:,'k), ^(a;,k) and i3(w,k) are all of the form 

D{uj,'k) = - iujDi{ui,'k) (41) 

where Dfi{uj,]i.) and -D/(a;,k) are are real-valued even functions of uj and k. The 
real parts Dji and Bfi must have a non-zero limit as a; —>■ 0 and k —> 0. All other 
parts of A, B and D must have finite limits as w —>■ 0 and k —?> 0. In the small w 
and k limit, the response function becomes 

, 1 \ ^ 5oo(0) , 

where we defined the diffusion constant D = 11^(0,0). 

The pole structure of G™ dictates that in the small uj and |k| limit, SpB{t,x) = 
S{pBit,'x.)) obeys the diffusion equation 

dtSpB = DW^5pB (43) 

This is our first example of a dissipative hydrodynamic equation. The conservation 
law, current algebra, thermodynamic stability, and the general analytic structure 
of the correlation functions are all the ingredients one needs to get this diffusion 
equation for baryon density. Hence diffusion is a very general phenomenon whenever 
there is a conserved current. Microscopic dynamics only enters through the value 
of the diffusion constant. 

If we now go back to the conservation equation 

dt^pB = -O^SJb (44) 

we can see that the diffusion equation above is equivalent to the constitutive rela¬ 
tionship 

6B = V^ = Dd^5pB (45) 

valid in the fluid cell rest frame. In the more general frame boosted by u^, it 
becomes 

= DAf^’^d^pB (46) 

where again A^’' = is the local 3-metric. 

The diffusion constant can be calculated by taking the appropriate limits of Gl 


lim lim —ImGL(w,k) = Zlgoo(O) 
w—>-0 k—>-0 OJ 


(47) 


which is our first example of the Kubo formula which relates an equilibrium corre¬ 
lation function to a dissipative coefficient. As the response function is singular in 
the small w and small k limit, the order of limits in the Kubo formula is important. 
The k —> 0 limit must be taken first. Since the transport coefficients are defined as 


Lorentz scalars, Eq.(47) is to be evaluated using the underlying microscopic theory 


such as thermal QCD in the rest frame of the equilibrium system. 
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3.3. Stress-Energy Tensor Correlation Functions 

To carry out the linear response analysis of the energy-momentum currents and 
their hydrodynamic and dissipative behavior, one need to know the Ward identity 
among the correlation functions. Defining the retarded correlation functions of 
TiJ.i' turned out to be not so straightforward due to the fact that the conserved 
quantities = f (PxT^^ are also the generators of the space-time evolution. 
Hence, in general the equal time commutators of are non-zero unlike the net 
baryon current case. 

To begin the analysis, consider the static partition function given by 


ZE[gE] = J 


-Se [0,ffE] 


(48) 


where 


SE[(l},gE] = J (Px J dr C{(j),gE) 


(49) 


is the Euclidean action and r is the imaginary time. Here, (j) denotes the collection 
of field variables and gE is the Euclidean metric. 

Using the Hilbert definition of the stress-energy tensor density, we have 


(r'^‘^(x))^ = -2 

The two point functions are given by 


^glu{x) 


\nZE[gE\ 


(50) 


{XE. Ve) = {TrT^^%XE)T^^{yE))E 

= 4 -, , - , , , -r In Ze[ gE] 


(51) 


^ ^gal}iyE)6g^AxE) 

where xe = {t, x) and %■ is the time ordering operator in t. Eor the tensor density 
, the covariant conservation law is 


d^{TnE + Kp{T^^)^ = Q 


(52) 


where T^^ = - g^p.p) is the Christoffel symbol. 


By differentiating Eq.(52) once more with respect to (/^^(j/), we obtain the Ward 
identity among the Euclidean correlation functions in the flat spac^^Ulll where g^ = 
5^^'' 




(53) 


Here k'^ = {lOn, k) and a;„ = 2TmT is the Matsubara frequency. 

To obtain the real time correlation functions, we perform the analytic continu¬ 


ation. From Eqs.(27) and (29) one can see that the analytically continuation 

uj^ —y — ik^ 6 


(54) 


changes the Euclidean correlation function to the real-time retarded correlation 
function. This also means that each time index of gets a factor of (—i). Since 
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our Minkowski metric is mostly negative, going from the Euclidean metric to the 
mostly negative Minkowski metric means 5p,u —>■ —g^^. The real-time version of 
Eq.(53) is then 

0 = ka 






(55) 


where = (a;,k). 


The presence of the single stress-energy tensor average terms in Eq. (55) implies 

is not the same as 


that the correlation function G' 


R 


'-'R 




eq 


(56) 


but differs by terms containing S{x — x') (as well as the contact terms containing 
spatial derivatives of 6{x — As the response function in the linear response 

theory, these delta-function terms do not matter since t and t' < t can never be the 
same. 

Defining the real-time correlation function by the analytic continuation of the 
Euclidean correlation function enables us to gain the following important relation¬ 
ship between the two 


G^"’“'^(0,k) = G^"’“^(0,k) 

which has a well-defined limit {p, and v here are not summed) 


^im G"(0,k) = ^im G"(0,k) > 0 


(57) 

(58) 


3.4. Momentum Diffusion and Shear Viscosity 

In this and the following section for the bulk viscosity, we will not consider finite 
net baryon density for simplicity. Eor analysis with finite ys, see Ref.^ 

Suppose that we set up a system where the the flow velocity at t = 0 has 
a single non-zero component in the a:-direction Ux{y) which depend only on y. 
In this situation, two layers of the fluid at y and at y -I- Ay have different fluid 
velocities in the orthogonal a;-direction. In ideal hydrodynamics, this difference is 
maintained because there is no dissipation. In a normal fluid, however, particle 
diffusion between the two layers will eventually make them move with the same 
equilibrated speed. How fast two layers equilibrate depends on the size of the 
scattering mean free path, which in turn determines the diffusion constant, or the 
shear viscosity. 

To set up a shear flow, let the perturbing Hamiltonian be 

6H{t) = -J dVe^*f"°(t,x)/?,(2/) (59) 

The corresponding linear response is 

/ OO 

dt'e{-ty^*'Gf’^°{t-t',ky) 

-OO 


(60) 
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for t > 0. 


In Eq.(551, appears in the following sequences when k = (0, ky, 0) 




xO,xy / 

R \ 

(w, ky) = ky) + P) 


^{G'^r (‘*^1 ky) + e) = fcyGfl ky 


Combined, these become 






(61) 

(62) 


(63) 


Except the extra P and e, the structure of these equations is exactly the same as 
the baryon current case. The following analysis is therefore a repeat of that case. 

In the w —>■ 0 limit, G^’^^{u!,ky) must have a well defined limit since it is a 
thermodynamic quantity. Furthermore, the imaginary part of G^’^^(uj,ky) must 
be an odd function of w. As in the baryon current case, the ky ^ 0 limit of the 
correlation functions must be well-behaved. Thus, we can parametrize as 

+ griky) + iujAxioj, ky)) 


GT^yiuj,ky) = 


k"^ — iw/DriuJ, ky) — uj'^Bt{oj, ky) 


and 




klie + griky) + iujAT{u},ky)) 


— iuj/DT(oJ, ky) — u}'^Bt{(jJ, ky) 


-P 


— s 


(64) 


(65) 


where gxiky) = G))°’^°(0, ky). Here the functions At, Bt and Dt all have the form 


Dt(oj, ky) = Dt{ui, ky) — iujDTioJ, ky) 


( 66 ) 


where D^{uj,ky) and D^p{uj,ky) are real-valued even functions of uj and ky. The 


.’rpyiu, rx.yj CUlU J^rpyM,n,yl 

real parts Dkp and Bkp must have a non-zero limit as w 


0 and —)■ 0. All other 


parts of At, Bt and Dt must have finite limits as w —>■ 0 and fcy —> 0. 

In the configuration space, the constant —e term becomes —e5'^{x — x'). In 
Eq.(60), this (5-function term does not contribute. Hence, in the small uj and ky 
limit, the evolution of is determined by iu: = DTky or 


{dt-DTdl)T-'^{t,y)=Q 


(67) 


where we defined the momentum diffusion constant Dt = D!p{0,0). This is our 
second dissipative hydrodynamic equation. The diffusion equation combined with 
the conservation law implies the constitutive relationship 

T^y(t, y) = DTdyp^^it, y) = gdyu^ (68) 

valid in the local rest frame. Here f] = DT{e + P) is the shear viscosity. It is clear 
from Eq.(67) that Dt has the physical interpretation of the diffusion constant for 


momentum diffusion. 
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Recognizing Eq.(68) to be a part of the spin 2 component of the second rank 
tensor, we can generalize this result to 

7rJ:^g(t,x) = ry (69) 

again in the rest frame of the fluid cell. In the moving frame, this becomes 

x) = = 2rycr'"‘' (70) 

where a^'' is the velocity shear tensor and is the the spin-2 projector defined 

by 


A^; = 1 (^AJiA^ + A:(A^ - 


(71) 


Here the label NS indicates that this is the Navier-Stokes form of the shear tensor. 
The Kubo formula for the shear viscosity is 


1 


lim lim —lTi\G^^'^^{(jj,ky) = Dt{s + P) = rj 


cl »—>-0 ky—^0 UJ 


(72) 


where we used the fact that (7 t(0) = P which can be determined from Eq.(55). 

The Kubo formula for the shear viscosity ry can be also expressed in terms of 
the full shear-tensor correlation functiorE^Hl 


11 = lim lim -ImGp*^ 

ccJ^O k^o lOw ^ 


'(w,k) 


where the shear-tensor is given by 




= T,, - {S.j/3)Ti 


(73) 

(74) 


3.5. Sound Propagation and Bulk Viscosity 


So far, only the diffusion type of hydrodynamic flow is discussed which are not the 
main bulk excitation. To get the main excitation which must also include the ideal 
hydrodynamics part, one needs to look at the disturbance in the energy density. 
This bulk excitation, of course, is the sound wave. 

Suppose we perturb the energy density with 


5H(t) = — J d^x (t, x)/3o (x) 

The linear response is then 

5(T00(t,k)) =/3o(k) 


dt'0{-t')e^*'G°J°°{t-t','k) 


Applying the conservation law to each index of G^’°‘^ in Eq.(55), we get 

k) = ui^e — a;^k^(e -|- P) -f k^GL(w, k) 

where 


(75) 


(76) 


(77) 


k^GUuj, k) = hkjkikm (g)(’'™(w, k) -t- (78) 
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In the small oj limit, Eq. gives 

Gl (w, k) « ^ (g™’™ (0, k) - (79) 

We also know that the imaginary part of Gi, must be an odd function of uj and the 
correlation functions are well-behaved in the k —^ 0 limit. The most general form 
consistent with these conditions is 


GL{uj,k) = 


(e + P- 


^Q{u},k)) 


k^ — ijp jZiuj, k) -I- iuj^R{uj, k) 
Here Z{uj,k), Q{uj,k) and i?(a;,k) all have the form 

Z{u}, k) = Zji{uj^ k) — iojZj{uj, k) 


(80) 


(81) 


where Zi^{uj,k) and Zj(u!,k) are real-valued even functions of w and k. The real 
parts Zfi and Rpi must have non-zero limits as w —> 0 and k —^ 0. All other parts 
of Z, Q and R must have finite limits as w —)■ 0 and k —)• 0. Matching the small w 
limit (79) demands that 


^ii(0,k) = 


e -I- P 


^ 00,00 ^ 


(82) 


(0,k) - e 

Up to the quadratic terms in w and k, the poles of Gl(w, k) for small w and |k| 
are determined by 


- Zr{0, 0)k2 -h iwZiiO, 0)k2 = 0 


(83) 


This has the structure of the dispersion relationship of a damped sound wave. 
Hence, in the small w and the small |k| limit, Zfl(0, 0) = is the speed of sound 
squared and ^/(O, 0) is the sound damping coefficient. 

To relate Zi{0, 0) to the shear and the bulk viscosities, let us consider the con¬ 
stitutive relationships once again. From the shear part, we already have the spin-2 
part of the stress tensor in the fluid cell rest frame 

diT‘°] (84) 


7rNs(^>x) = Dt 


QirpjO gjj^iO _ ‘^9 ^ 

3 


To this we add a spin-0 part to get 

x) = Dt + jg^^diT‘° (85) 

The energy conservation law in the local rest fram^now becomes using Eq.(14) 
with the dissipative part given by Eq.(85) 


0 = d^d.T^'' 

= die - - Dr^'^^dte - 

(-OJ^ + ^^^k^ - *(4 U>t/ 3 -b i)k^uj) Se (86) 

^ In the local rest frame, u(t, x) = 0, but diUj 7 ^ 0. 
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where we used g^^didj = — and dte 
identify = dP/de = and 


—diT'‘°. Comparing with Eq.(83), one can 


Zj(o,o) = r = 4i:)T/3 + 7 


(87) 


as the sound attenuation constant. Since we have already identified Dt = r]/(e+P), 
this allows us to identify 7 = (^/{e + P) where C is the bulk viscosity. In the 
fluid cell rest frame, the added term corresponds to the constitutive relationship 
IIns = —(diU^. In the general frame, this becomes 


IIns = 


( 88 ) 


Again, the label NS indicates that this is the Navier-Stokes form of the bulk pressure. 
The minus sign in Eq. ( 88 ) makes sense since the effective pressure P +11 should be 


less than the equilibrium pressure when the fluid is expanding (positive diU^). 

We have so far identified .^^^(0,0) and Zj{0,0) as the speed of sound squared 
and the sound attenuation coefficient. The role of R{uJ, k) is still to be identified. 
In the Kubo formula for the attenuation coefficient, i?/j(0,0) appears as 

Gi(a;,k) 


lim lim Im ■ 

cj —>-0 k —>-0 


= (e + P){Zii0, 0) - Zfl(0,0)2 Rr{0, 0)) 


(89) 


which does not allow one to identify the right hand side with T if 77^(0, 0 ) ^ 0 . 
Actually, the left hand side of Eq.(89l does yield T. It is just that we have not 


been consistent in power counting since the wave equation Eq.( 86 ) does not contain 
0{uj^) term while the pole of does. One may consider this discrepancy as the 


first sign of the trouble with the first order constitutive relationship Eq.(85). 


3.6. Second Order Viscous Hydrodynamics 

Let us consider the consequence of having the first order constitutive relationship 


more closely. The diffusion equation with a source S 


{dt — DtV^) 6n{t, x) = S{t, x) 

(90) 

has the solution 


Sn{x) = j d'^x' Gpi{x — x') S{x') 

(91) 

Here the retarded Green function is 


p AD{t-t') 

Gh(x x') = «(i 

(92) 


If one has a point source, S{x') = NoS{x'), then 6n{t,x) = NoGii{t,x). At t — 0, 
the space is empty except at the origin. But at any time after that, there is non-zero 
Sn everywhere. This is clearly acausal. 

On the other hand, the solution of the sound equation 

- {dt - = Six) 


(93) 
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for a point source S{x) = Ao(5(a;) is 

Se{t, x) = AoO{t) 


1 (5(|x| - Vst) 
47 r |x| 


(94) 


This is causal since the disturbance only moves with the speed of sound. 

The origin of acausality in diffusion is the mismatch between the number of 
time derivatives and the number of spatial derivatives in the diffusion equation. 
The diffusive dispersion relationship uj = —iDh^ gives the group velocity 




(95) 


which becomes large in the large k limit. This problem can be remedied if one 
replaces the constitutive equation, J* = Dd’'n with a relaxation type equation 


dtr = - — {.r-Dd^n) 
tr 


then the conservation law becomes 

—dtn^ — 

Tr th 

For large k where we previously had a problem, we now have 


d?n = —didtJ^ = ——dtn + — 




2 vie 


(96) 

(97) 

(98) 


with the propagation speed vr = y/Djrl. 

This type of relaxation equation was actually anticipated already: Up to the 


second order in w, the poles of the the density-density correlator in Eq.(40) are 
determined by 

D\e — iuj — eDB = 0 


Comparing, we see that 


B = tr/D 


(99) 


( 100 ) 


For the viscous stress-energy tensor components, the following relaxation equations 
apply in the local rest frame 


-I- — ^ TT®-^ = — TT®-^ 


— j - —'‘NS 

dt -\ -) n = —IIns 

TnJ Tn 


where 


^Lm rjiLm 

IT = 1 — 


-Tt 


is the traceless part of the stress tensor and 


1 , 


n = - -T, 


p 


( 101 ) 

( 102 ) 

(103) 

(104) 
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is the bulk pressure. They are not to be identified with the Navier-Stokes forms 


(84) and ( 88 ). Rather, they will relax to the Navier-Stokes forms. 


To see if the acausality in the momentum diffusion is cured, we start with the 
momentum conservation 


= -diT 


M 


Applying the curl gives 


dtT^T = -^ijkdjdiTT 


kl 


(105) 


(106) 


where we defined = CijkdjT^^. Applying {dt + 1/t^) and using Eq.(84) yields 

0 = + dt- Dt^^) (107) 

As long as < 1, this is now causal. 

For the sound modes, we start with the conservation law in the local rest frame 


a?e = -b - v^n 


(108) 


Applying (r^i9t-|-l)(Tn9t-|-l) to Eq.(108) and using Eqs.(lOl) and (102), one obtains 


the following dispersion relation for the bulk mode propagation (in this case 6e) 

n 4 2 2i 2 ADt ,22 i 2 2 

0 = TttTuUJ — TTrTnVgUJ K — Tn—-—k W — Tjryk w 

2 I 2i 2 
—UJ -|-4^>,k —t 


“■ +7 + v'^{t^ + rn)^ k^w -b iu}^ + Tn) 


(109) 


Comparing with the small uj and small |k| expansion of the denominator in Eq.(80), 
we can identify Zfi{0, 0 ) = 


Z/(0, 0) = + 7 + vUt 77 + m) 


and 


-Ri?(0,0) = (r,, -bTn)/u^ 

The Kubo formula for the damping constant now makes more sense 
GL(a;,k) _ ^ ^ 2 


( 110 ) 

( 111 ) 


lim lim Im - 

uj—>-0 k—>-0 CJ 


= (e + F)(Zj(0, 0) - Zn(0, 0)" Er(0, 0)) 


= (e + P) f + 7 


V 3 




and for the bulk viscosity only. 


( = lim lim — flmG'i(w,k) — -Im Go^’'^^(a;, fc„) 


( 112 ) 


(113) 
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The Kubo formula for the bulk viscosity ( is also available in terms of the pressure- 
pressure correlation functiorUnHIl 


pp, 


C = lim lim — ImGp 
o;—>-0 k—vO LU 


w,k) 


(114) 


Using the fact that the correlation functions of vanishes in the k —0 limit, 
one can use in place of P the trace or the combination P — v^e to make it 

more explicit that the bulk viscosity is non-zero only if the conformal symmetry is 
broken. The Kubo formulas for the relaxation times and rn are not simple to 
determine in this analysis. Simple Kubo formulas for has been worked out in 
Refs Emu as 


= — lim lim 

uj^O k—>-0 2 

although a simple Kubo formula for Tn is still to be found. 


(115) 


In the dispersion relation Eq.(109), the 4-th order terms of 0{uj‘^) and 0(a;^k^) 
are present but O(k^) terms are not. This may seem unsatisfactory since there is no 
reason why this term should be small compared to the other two 4-th order terms 
when Lo ^ Us|k|. However, one should recall that hydrodynamics is valid only in the 
long wavelength and small frequency limits. From this point of view, the 4-th order 
terms are not so important. They become, however, significant when the equations 
are solved numerically that can include short wavelength excitations. Fortunately, 
this Israel-Stewart form of second order hydrodynamic^^ (comprising of Fqs.(96l, 
(1011 and (102)) is shown to be stable in RefsElEll 

If one wants to include 0(k‘^) terms, then one can modify the relaxation equa¬ 


tions (101) and (102) to include second derivatives of . However, doing so not 


only generates O(k^) terms, but it also generates (incomplete) terms involving 5 
and 6 factors of uj and k. Since higher order terms begin to matter at large w and 
|k|, having higher and higher order of frequency and momentum (or derivatives in 
the configuration space) does not guarantee that the numerical solution in this limit 
becomes more and more faithful to the real spectrum. One just needs to be careful 
not to interpret high frequency and momentum modes as physical. 

As for the calculation of the viscosities, full leading order perturbative QCD 
results for both the shear viscosity and the bulk viscosity have been obtained in 
Refs.iSnHH] 

using the Kubo formulas illustrated above. QCD is an asymptotically 
free theory.l^^H^ In principle, it admits a perturbative expansion only when the 
energy scale exceeds at least a few GeV. Since the typical QGP energy scale is 
less than 1 GeV, the strong coupling is not small. Phenomenologically, we must 
have as « O.SEIEZI This value may look weak, but the gauge coupling itself g = 
y/Airas « 2 is not so small. In the perturbative many-body QGD, g (or gj^'P) is 
the expansion parameter not Therefore although the analysis performed 

in Refs.innHH] 

are nothing short of tour de force, having g k 2 makes numerical 
values obtained in perturbation theory not too reliable. At this point, reliable 
first principle calculations at large g can only be performed on numerical lattice 
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in Euclidean space. Lattice QCD can straightforwardly compute static properties 
such as the equation of state. However, calculations of dynamic properties such as 
the viscosities become more complicated as they involve estimating real continuous 
functions from a finite set of discrete Euclidean data. Nonetheless, a great deal 
has been accomplished in computing the properties of QGP through lattice QCD 
calculationd^^MSl gg through effective models such as the hadron resonance 

gas model (HR)HI1^ and the AdS/CET correspondenceH^HlEl 

The purpose of this section has been to show that the hydrodynamics is very 
general. No matter what the system is, there usually is a regime where hydro¬ 
dynamics is in some way applicable as long as there exist “macroscopically small 
but microscopically large” length and time scales. For more detailed analysis of 
the length and time scales and also for demonstrating more general structure of 
hydrodynamic equations, we now turn to the kinetic theory. 


4. Hydrodynamics from kinetic theory 


4.1. Length scales and validity of hydrodynamic approximations 

Kinetic theory describes a medium microscopically, by following the evolution of 
the phase-space distribution function f{x^p), a Lorentz scalar that describes the 
probability of finding a particle with four-momentum at space-time position 
x". Classical kinetic theory assumes that the particle momenta are on-shell, p^ = 
TO^, which requires the system to be sufficiently dilute and the mean free paths 
sufficiently long to ignore collisional broadening effects on the spectral function 
p{p) = 2'K5(jP‘—m?') that defines the particles’ propagator. The defining equation 
of classical kinetic theory is the Boltzmann equation. 


P^dfj,f{x,p) = C{x,p), 


(116) 


where C{x,p) is the collision term in which the strength of the interaction enters 
through their scattering cross sections. Especially for massless degrees of freedom, 
its detailed form can be quite complicated.^^ A popular simplification of the collision 
term is the relaxation time approximation (RTA)|^ 


C{x,p) 


P^u^{x) 

Tre\{x,p) 


feq{x,p)-f{x,p) 


(117) 


where the relaxation time Trei in general depends on position through the local 
density and can also depend on the local rest frame energy of the particles (indicated 
by the p-dependence). Classical kinetic theory is valid if this relaxation time Trei, 
and the associated mean free path Amfp = {ip/E)xiei), are sufficiently large. In 
other words, interactions among the constituents must be weak. 

‘^The Boltzmann equation with this RTA-approximated collision term is known as the Anderson- 
Witting equation. 
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Hydrodynamics is valid if the system is close enough to thermal equilibrium that 
its local momentum distribution (and therefore its macroscopic fields, such as par¬ 
ticle and energy density and pressure, which can all be expressed as moments of the 
local momentum distribution) can be characterized by a small number of thermo¬ 
dynamic and transport parameters, such as temperature, chemical potential, shear 
and bulk viscosity, etc. This requires efficient interactions among the constituents 
of the medium because otherwise any kind of macroscopic dynamics involving lo¬ 
cal expansion or compression or shear of the fluid will drive its local momentum 
distribution away from its near-equilibrium form. Hydrodynamics works best for 
systems made of strongly interacting constituents. 

Does this mean that the validity of kinetic theory and hydrodynamics are mu¬ 
tually exclusive? Not necessarily. To gain clarity consider a relativistic system of 
(almost) massless degrees of freedom. It can be characterized by three length scales, 
two microscopic and one macroscopic one: 

• the thermal wavelength Ath ~ 1 /T 

• the mean free path Amfp ~ l/((c’''i')^) where (av) is the momentum-averaged 
transport cross section times the relative speed (« 1 in units of c) of the 
colliding objects, and n is the density of scatterers 

• the length scale Lhydro over which macroscopic fluid dynamical variables 

vary; it can be defined in many ways that give quantitatively different but 
similar order of magnitude results: ~ ~ 

The ratio between the two microscopic scales characterizes the magnitude of the 
transport coefficients rj (shear viscosity), ( (bulk viscosity), and k (heat conductiv¬ 
ity): 

-^mfp ^ 1 1 ^ 1 1 (118) 

Ath {f^)n Ath (o-)Ath s s’ s’ s’ 

where we used rj, C, Tn ~ l/((CT)Ath) ~ AmfpT'’^ and the entropy density s ~ 4n ~ 
for a near-thermalized system of particle number density n for massless degrees 
of freedom. 

In terms of the two microscopic length scales we can define three regimes of 
microscopic dynamics: 

(1) Dilute gas regime: 

(119) 

Ath s 

This is the weak-coupling regime where the microscopic system dynamics can be 
described in terms of on-shell quasi-particles and many-body correlations are 
suppressed. In this regime the Boltzmann equation applies. 

(2) Dense gas regime: 

^ (120) 

Ath S 
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In this case interactions happen on the scale Ath- We call this the moderate cou¬ 
pling regime where the microscopic system dynamics must be described by off- 
shell quasiparticles (whose spectral functions have a finite collisional width) and 
many-body correlation effects are non-negligible. Here the Boltzmann equation 
must be replaced by a quantum kinetic approach based on Wigner distributions, 
and the BBGKY hierarchy of coupled equations for the iV-body distribution 
functions can no longer be efficiently truncated. 

(3) Liquid regime: 

^--«1 ^ (121) 
Ath s 

This is the strong-coupling regime where the system has no well-defined quasi¬ 
particles and no valid kinetic theory description. 

To judge the validity of a macroscopic hydrodynamic approach we compare the 
microscopic to the macroscopic length scales. To simplify the discussion, let us 
agree on using the inverse of the scalar expansion rate 6 — to represent the 

macroscopic length scale Lhydroj^ The figure of merit controlling the validity of a 
fluid dynamic picture is the Knudsen number. 

Kn = A„,fp • e ~ ^Ath Ot,,,. (122) 

The Knudsen number is the small parameter that controls the convergence of the 
expansion in gradients of thermodynamic quantities that underlies the derivation 
of hydrodynamics as an effective theory for the long-distance dynamics of a general 
quantum field theory.^ Again, we can use it to define three regimes: 

(1) Ideal fluid dynamics: 

Kn « 0 « 0 or 6* « 0 such that Or^^i ~ 0 (123) 

s 

(2) Viscous fluid dynamics: 

Kn < 1 <;=^ - 01 9 small such that OT^ei ^ 1 (124) 

(3) Hydrodynamics breaks down: 

Kn 1 <;=^ - or 0 large such that 0Trei ^ 1 (125) 

s 

In high-energy heavy-ion collisions, the initial energy deposition occurs in an 
approximately boost-invariant fashion along the beam direction, leading to an ex¬ 
pansion rate 9 that diverges like l/r for small time r after impact]^ On the other 

®The shorter L^ydroj the faster the system is driven away from local equilibrium. The scalar 
expansion rate directly drives the bulk viscous pressure 11. It is par ametrically of the same order 
as the shear tensor defined in Eq. |70l that drives the shear viscous 

pressure and as the diffusion force /^ = V^(/l/T) associated with space-time gradients of 
conserved charge densities that drives the heat flow 

tflere, r = is the longitudinal proper time and the boost-invarian ce re fers to independence 

of the space-time rapidity r] = tanh“^( 2 :/f). For more details, see section 5.1 
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hand we now know that the quantum mechanical uncertainty relation places a lower 
bound rj/s > l/(47r) on any system, even at infinitely strong coupling EESUl There¬ 
fore, hydrodynamics is inapplicable during the earliest stage of a heavy-ion collision. 
At the end of a heavy-ion collision, the mean free paths of hadrons become large 
compared to the Hubble radius ^ 1/0 oi the expanding fireball, and hydrodynamics 
breaks down again. This process is called kinetic decoupling. Between the early 
pre-equilibrium and the final decoupling stage stretches an extended period of ap¬ 
plicability of viscous fluid dynamics. The most important factor ensuring this is 
the strong collective coupling of the quark-gluon plasma (QGP) phase which is 
characterized by a small specific shear viscosity rj/s ~ ( 2 — 3 )/( 47 r).ISSHlll 

Note that the validity of hydrodynamics does not rely directly on 77/s being 
small, only on (rj/s) • {O/T) being small. So, strictly speaking, strong coupling is 
not required for hydrodynamics to be valid. Only in extreme situations, such as 
heavy-ion collisions which are characterized by extreme expansion rates, does hy¬ 
drodynamics require very strong coupling. In this case, hydrodynamics is applicable 
even though classical kinetic theory is not, because very strongly coupled quantum 
field theories do not allow a description in terms of on-shell quasi-particles. It is gen¬ 
erally believed that the very earliest stage of a heavy-ion collision has no well-defined 
quasiparticles at all and is better described by a theory of classical or quantum fields 
than by a (quantum) kinetic approach. On the other hand, weakly coupled system, 
with very large values of 77/s in which the applicability of (even classical) kinetic 
theory is ensured, can still be describable macroscopically through fluid dynamics 
if they are sufficiently homogeneous and expand slowly. In this case the smallness 
of 6 can compensate for the largeness of 77/s, resulting in a small Knudsen number. 
Systems with a large 77/s but a small product (77/s) • (9/T) admit simultaneous 
microscopic classical kinetic and macroscopic hydrodynamic descriptions. In the 
following subsections we will study such systems to derive the macroscopic hydro- 
dynamic equations from the microscopic kinetic theory. Hydrodynamics being an 
effective long-distance theory, the form of the resulting equations does not rely on 
the validity of the underlying kinetic theory (although the values for the transport 
coefficients do); they can therefore also be applied to a strongly-coupled liquid such 
as the QGP. 

4.2. Ideal fluid dynamics 

We define 77-moments of the distribution function weighted with some momentum 
observable 0{p) by 

{0{p)) = J 0{p)f{x,p) = g J j^^^^Oip) f{x,p) (126) 

{g is a degeneracy factor) and p° = Ep = \JrrE- -|- p^. The particle number current 
and energy momentum tensor are then written as 

j'^ = (77^), T^'' = lp^p''). 


(127) 
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Usually there is more than one particle species in the system, and the conserved 
baryon charge current Jg and energy-momentum tensor are given in terms of 
linear combinations of {p^)i and {p^p^)i where the subscript i labels the particle 
species whose distribution function is fi{x,p): 

= = E( 128 ) 

i i i i 

here bi is the baryon charge carried by each particle of species i. For simplicity, we 
restrict the following discussion to a single particle species. 

The particle number current and energy-momentum tensor take their ideal fluid 
dynamical form 


■U Li 

Jid = 






where the spatial projector in the local rest frame (LRF) 
we assume that the system is locally momentum isotropic: 

'p^u^[x) - p{x) 


f{x,p) = fiso{x,p) = fi, 

The local equilibrium distribution 

/eq(C) = - 


Tix) 


+ I 


(129) 

is given in Eq. 

(130) 

(131) 


where ( = {p^Ui_i{x)—p{x))/T{x) and a = 1,—1,0 for Fermi-Dirac, Bose-Einstein, 
and classical Boltzmann statistics, respectively, is a special form of fiso{x,p). It is 
defined as the distribution for which the collision term C(x,p) in the Boltzmann 
equation (116) vanishes. Note that the ideal fluid decomposition ( |129 ) does not re¬ 
quire chemical equilibrium, i.e. it holds for arbitrary values of the chemical potential 
p{x), nor does it require complete thermal equilibrium, i.e. /iso is not required to 
depend on its argument exponentially as is the case for the equilibrium distribution 
(131). If the dependence is non-exponential, the collision term in the Boltzmann 


equation is non-zero, but its p^-moment still vanishes, J^p'^C = 0, due to energy- 
momentum conservation. 

The ideal hydrodynamic equations follow by inserting the ideal fluid decompo¬ 
sition (129) into the conservation laws Eq. (0: 


h = —nO, £ = —{£+P)9, 


\jt^p 


e+P 1 + ci 




(132) 


where the very last expression assumes an EOS of type P = c^e. F denotes the LRF 
time derivative of a function P, P= Pt-P = m^ 9^P, and the spatial 

gradient in the LRF. Thus, 9^ = -I- 

Equations (132) can be solved numerically for the local particle density n{x), 
energy density £{x), and flow velocity u^{x), with the temperature T{x), chemical 


SNote that in curvilinear coordinates or curved space-times, the partial derivative must be 
replaced by the covariant derivative 
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potential fj,{x) and pressure P{x) following from the equation of state (EOS) of 
the fluid. Local deviations from chemical equilibrium result in a non-equilibrium 
value of the local chemical potential fJ-{x) and a non-zero right hand side in the 
current conservation equation for j^. Deviations from thermal equilibrium (while 
preserving local isotropy) must be accounted for by a non-equilibrium pressure in 
the EOS P{s,n). In both cases, the conservation laws, Eqs. Q, lead to a non¬ 
vanishing entropy production rate 9^5'^ ~ l/TVei 


4.3. Viscous fluid dynamics 

4.3.1. Navier-Stokes (NS), Israel-Stewart (IS) and Denicol-Niemi-Molnar- 
Rischke (DNMR) theory (VHydro^ 


Israel-Stewart (IS) and Denicol-Niemi-Molnar-Rischke (DNMR) second-order vis¬ 
cous fluid dynamic^^ are obtained by using in (127) for f{x,p) the ansatz 


f{x,p) = /is 


p^Uf^{x) - y{x) 
T{x) 


6f{x,p). 


(133) 


The correction Sf describes the deviation of the solution f{x,p) of the Boltzmann 
equation from local momentum isotropy. It is supposed to be “small”, in a sense 
that will become clearer below, and will thus be treated perturbatively. 

Most authors set /iso = feq, i-e. they expand around a local equilibrium state. 
To obtain the correct form of the hydrodynamic equations this is not necessary; 
only the form of the equation of state P{e,n) and the values of the transport 
coefficients depend on this choice. We, too, will make this choice for simplicity, but 
emphasize that under certain conditions the perturbative treatment of Sf may be 
better justified if the leading-order distribution /iso is not assumed to be thermal. 

For later convenience we decompose p^ into its temporal and spatial components 
in the LRF: 


(134) 

where Ep = Up,p^ and p^^'^ = A^'^pi, are the energy and spatial momentum compo¬ 
nents in the LRF. Then 


= {Ep), e={El). 


(135) 


The decomposition (133) is made unique by Landau matching: First, define 
the LRF by solving the eigenvalue equation ([^ with the constraint u^Up = 1 which 
selects among the four eigenvectors of the timelike one. Eq. 0 fixes the flow 
vector u^(x) and the LRF energy density. Next, we fix T{x) and pl(x) by demanding 
that Sf gives no contribution to the local energy and baryon density: 


{Ep)s = {El)s = 0 . 


(136) 


Inserting (133) into (127) we find the general decomposition 




(137) 
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with a non-zero number flow in the LRF, 

= (138) 

a bulk viscous pressure 

n= - (139) 

and a shear stress 

= (140) 

where {■ ■ ■)s indicates moments taken with the deviation 6f from /iso- In the last 
equation we introduced the notation 




(141) 


where is the spin-2 projector introduced in Eq.(71) denoting the traceless and 
transverse (to u^) part of a tensor . The shear stress tensor thus has 

5 independent components while which is also orthogonal to by construction, 
has 3 independent components. 

Using the viscous hydrodynamic decomposition (1371 in the conservation laws 
= 0 and = 0, we obtain the vHydro viscous hydrodynamic evolution 
equations 


n = -n9 - 
i = — (£-\-P 4-11)0 -f 

ie+P+U)uf^ = Vf^iP+U) - -4 Tr>^''u^. (142) 


where is the velocity shear tensor introduced in Eq.(70). They differ 

from the ideal fluid dynamical equations (132) by additional source terms arising 
from the dissipative flows. Altogether, the deviation 6f has introduced (3-|-l-|-5)=9 
additional dissipative flow degrees of freedom for which additional evolution equa¬ 
tions are needed. These cannot be obtained from the macroscopic conservation 
laws but require input from the microscopic dynamics. In a system that is initially 
in local equilibrium, the deviation 6f is caused by the dynamical response of the 
system to gradients in the thermodynamic and flow variables. The forces that drive 
this deviation can be classified by their Lorentz structure as a scalar, a vector and 
a tensor force: 


scalar force: 6 = (scalar expansion rate); 

vector force: (fugacity gradient); 

symmetric tensor force: (velocity shear tensor), 

antisymmetric tensor force: V'^u^) (vorticity tensor). (143) 
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These forces generate dissipative flows, the scalar bulk viscous pressure 11, the heat 
flow vector V^, and the shear stress tensor The strength of the forces (143) 

driving the system away from local equilibrium is characterized by the Knudsen 
number. The system response can be characterized by inverse Reynolds numbers 
associated with the dissipative flows® 


|n| 

— p > 


Ry — 




p 


R-^ = 


V' 




P 


(144) 


Due to the time delay r^ei between the action of the force and the system response, 
built into the collision term of the Boltzmann equation, the inverse Reynolds num¬ 
bers are not necessarily of the same order as the Knudsen number: For example, 
an initially small bulk viscous pressure can remain small, due to critical slowing 
down, as the system passes through a phase transition, even though the bulk vis¬ 
cosity becomes large during the transition.® Conversely, strong deviations from 
local equilibrium during a rapidly expanding pre-equilibrium stage in heavy-ion 
collisions can lead to large initial values for the dissipative flows, and a slow equili¬ 
bration rate may cause them to to stay large for a while even though the expansion 
rate decreases with longitudinal proper time as 1/r. Deviations from equilibrium, 
and the accuracy of their description by viscous fluid dynamics, are therefore con¬ 
trolled by a combination of Knudsen and inverse Reynolds numbers.® 

The 9 equations of motion describing the relaxation of the 9 dissipative flow 
components are controlled by microscopic physics, encoded in the collision term on 
the right hand side of the Boltzmann equation, and can be derived from approximate 
solutions of that equation. This was first done almost 50 years ago by Israel and 
Stewart in Ref.,® but when the problem was recently revisited it was founcGIEMlIlll 
that the relaxation equations take a much more general form than originally derived. 
Specifically, Denied et aZ.® found the following general structure 


Tnld -|- n — — (^0 -\- Sf -f ZC -f 7^, 

+J^^+ICf^+ 7 ^'^, 

-b -b (145) 


Here all calligraphic terms are of second order in combined powers of the Knudsen 
and inverse Reynolds numbers. terms contain products of factors that are each of 
first order in the Knudsen and inverse Reynolds numbers; 1C terms are second order 
in Knudsen number, and R terms are second order in inverse Reynolds numbers. 

In relaxation time approximation, with an energy-independent relaxation time 
Tj-ei, the relaxation times for the dissipative flows all agree with each other: 
Tn = Ty = T,r = TreiThe Same does not hold for more general forms of the col¬ 
lision term. If we set in Eqs. (145) the relaxation times and all other second order 


*'The energy-momentum tensor is symmetric, so dissipative flows have no antisymmetric tensor 
contribution, but the antisymmetric vorticity tensor couples to the other dissipative forces and 
flows at second order in the Knudsen and inverse Reynolds numbers. 
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terms to zero, we obtain the equations of relativistic Navier-Stokes theory: 

n = -C6», = kI^^, = 2r]a>^^. (146) 


The relaxation equations (145) have solutions that, for sufficiently small expansion 


rates (see below), approach asymptotically (at times r :g> Tn^yjr) the Navier-Stokes 


values (146). However, plugging the Navier-Stokes solutions (146) directly back 


into the decompositions (137) and using them in the conservation laws ([^ leads to 
viscous hydrodynamic equations of motion that are acausal and numerically unsta- 
bleElEll The physical reason for this is the instantaneous response of the dissipa¬ 


tive flows to the dissipative forces encoded in Eqs. (146) which violates causality. 


A causal and numerically stable implementation of viscous fluid dynamics must 
account for the time delay between cause and effect of dissipative phenomena and 
therefore be by necessity of second order in Knudsen and Reynolds numbers. 

The first relativistic causal second-order theory of viscous fluid dynamics was 
Israel-Stewart (IS) theory.l^ It amounts to dropping the K, and TZ terms in (145) 
and replacing (for massless particles) J —>■ — Irn^H, —>■ —TyOV^, and —>■ 

— . The importance of keeping specifically these second-order jZ-terms for 

the preservation of conformal invariance in a system of massless degrees of freedom 
was stressed by Baier et 

For conformal systems the resulting Israel-Stewart relaxation equations can be 
written in the fornPH 

11 = - 


Tv 






(147) 


with effective transport coefficients and relaxation times that are modified by the 
scalar expansion rate as follows: 

C , K, 


C' = 


l+7n ’ 


l+7v’ 


7?' = 


V 


1+771- 


l+7i 


{i = IV,V,TT), (148) 


where 7 ^ = for i = H, tt and 7 ^ = Oti for i = V. 

These relaxation equations describe an asymptotic approach of the dissipative 
flows to effective Navier-Stokes values that, for a positive scalar expansion rate 9, are 


reduced relative to their first-order values (146) by a factor l+ 7 i, while the effective 


rate of approach to this effective Navier-Stokes limit is sped up by the same factor. 
7 i involves the product of the scalar expansion rate and the respective relaxation 
time. So, compared to a more slowly expanding system, a rapidly expanding system 
with the same microscopic scattering cross sections is characterized by lower effective 
viscosities and shorter effective relaxation times 1^51^ 

When all the second order terms are kept, the DNMR equations become quite 
complicated. The J', K, and TZ terms listed by DNMRp^ add up to 16 terms for 
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n, 18 terms for V^, and 19 terms for each with its own transport coefficient. 
While many of these transport coefficients have been computed for massless the¬ 
ories at both weak and strong coupling (a subject too rich to be fully reviewed 
here, see papers by Denicol et al. jMHIlII Moore et and Noronha et al^ for 

lists of references to relevant work), quite a few are still unknown. Furthermore, 
the QGP is neither weakly nor strongly enough coupled, nor is it sufficiently con¬ 
formally symmetric for any of these calculations to be quantitatively reliable for 
heavy-ion collisions. For these reasons, many practical applications of viscous fluid 
dynamics employ phenomenological values for the transport coefficients, and work 
studying which terms need to be kept and which might be of lesser importance is 
still ongoing .1^ 

One important non-linear coupling mechanism that enters at second order are 
bulk-shear couplings where shear stress drives a bulk viscous pressure and vice 
versaUnMElIZll Heavy-ion collisions are characterized by initially very large differ¬ 
ences between the longitudinal and transverse expansion rates that cause large shear 
stress. The latter, in turn, creates a bulk viscous pressure via bulk-shear coupling 
that can dominate over the one generated a la Navier-Stokes by the scalar expansion 
rat^^ and may actually be able to flip its sign. This should be taken into account 
in phonemenological applications of viscous hydrodynamics to heavy-ion collisions. 


4.3.2. Anisotropic hydrodynamics (aHydroJ 

The dissipative flows are given by moments of the deviation Sf{x,p) of the distribu¬ 
tion function from local equilibrium, and their relaxation equations are derived from 
the Boltzmann equation using approximations that, in one way or another, assume 
that Sf is small. However, systems featuring strongly anisotropic expansion, such 
as the early evolution stage of the fireballs created in ultra-relativistic heavy-ion 
collisions, generate strong local momentum anisotropies: the width of the LRF mo¬ 
mentum distribution along a certain direction is inversely proportional to the local 
expansion rate in that direction, and this momentum-space distortion growth with 
the magnitude of the shear viscosity. In viscous hydrodynamics, where we expand 


/ around a locally isotropic LO distribution (see Eq. (133)), this local momentum 


anisotropy must be absorbed entirely by i5/, making Sf large and rendering the 
approximations used for calculating the evolution of the dissipative flow generated 
by Sf unreliable. Indeed, even for moderate specific shear viscosities rj/s ^ 5—10 
the (negative) longitudinal component of the viscous shear pressure can become so 
large in Israel-Stewart theory that it overwhelms the thermal pressure, resulting in 
a negative total pressure along the beam direction - which, according to the kinetic 
definition Pl = \{pI), should never happen. 

Anisotropic hydrodynamic^^3Zll is based on the idea to account already in the 
LO distribution for the local momentum anisotropy resulting from anisotropic ex- 
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pansion, by parametrizing^ 


f{.X,p) = fRs{x,p) = /iso 


/ ^/pf,^^^•'{x)p^ 

[ m 



(149) 


where = u^{x)u''{x) + ^{x)z^^{x)z'^{x), z^{x) being a unit vector in longitu¬ 

dinal 0 direction in the LRF. The subscript RS refers to Romatschke and Strickland 
who are the authors of Ref.^ This distribution is characterized by 3 flow parame¬ 
ters u^{x) and three “thermodynamic” parameters: the “transverse temperature” 
A(x), the effective chemical potential fj-{x), and the momentum-anisotropy param¬ 
eter /(a;). Inserting (149) into (127) yields the aHydro decomposition 


j^s = T^s = + [Pl - Pt)z^^z’', (150) 

nRs = {E)rs = P-o(C) »T-iso(A, /), ers = (P^)rs = P(0 eiso(A, jl), (151) 
Pt,L = (Pt,l)rS = 77T,L(0^iso(A,/i). 


For massless systems, the local momentum anisotropy effects factor out via the 
77(^)-functions!^ 




1 

VttI’ 

3 / l + (g^-l)p(0 \ 

2/1 / + ! 


= 2 

-- 


1 arctan 

3 / (/+l)P(/)-l 

/ V / + 1 


(152) 


The isotropic pressure is obtained from a locally isotropic equation of state 
Piso(A,/i) = Piso(£iso(A,/i), niso(A,/I)). For massless noninteracting partons, 
Piso (A, p) = |eiso(A,/i) independent of chemical composition. To compare with 
ideal and IS viscous hydrodynamics, we need to assign the locally anisotropic sys¬ 
tem an appropriate temperature T{x) =T(^^{x), A{x), fl{x)) and chemical potential 
fj,{x) = p{^{x), A{x), fl{x)), thinking of /rs(/, A) as an expansion around the locally 
isotropic distribution fiso(T). For this we impose the generalized Landau matching 
conditions £rs(/, A,/I) = £iso(T, Z^) and For example, using an exponential (Boltz¬ 
mann) function for /iso with p = p = 0, one finds T — AP^A(^)_ With this matching 
we can write 


Trs = + nRs)A'^'^ + 

AP -h Brs = -\ J PaA“'^Pi3(/Rs - /iso) (= 0 for m = 0), 

M. f ^ ^ ,x^x^ -2z>^z- 

= J P P T/rS-Ziso) = {Pt-Pl) —-—g-• 


(153) 

(154) 

(155) 


We see that 7r{(g has only one independent component, Pt—Pl, so aHydro 
leaves 4 of the 5 components of unaccounted for. For massless particles 
we have (Pt—PL)/P iso(e) =P-t(/)—Pl(/), so the equation of motion for 7r^g 
can be replaced by one for /. For (2-|-l)-dimensional expansion with longitudi¬ 
nal boost-invariance these equations can be found and were solved numerically 
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by Martinez et For m ^ 0 we need an additional “anisotropic EOS” for 

(AP/Piso) = {2Pt+Pl) / {3Piso) — 1, in order to separate AP from the viscous bulk 
pressure If. 


4.3.3. Viscous anisotropic hydrodynamics (VaHydro^ 

As explained above, AHYDRcdEH accounts only for one (albeit largest) of the five 
independent components of the shear stress tensor . It can therefore not be 
used to compute the viscous suppression of elliptic flow which is sensitive to e.g. 
^xx_^yy ^ On the other hand, since the four remaining components of the shear 
stress tensor never become as large as the longitudinal/transverse pressure differ¬ 
ence (with smooth initial density profiles they start out as zero, and with fluctuating 
initial conditions they are initially small), they can be treated “perturbatively” a 
la Israel and Stewart, without running into problems even at early times. Com¬ 
bining the non-perturbative dynamics of Pl—Pt via aHydro with a perturbative 
treatment of the remaining viscous stress terms a la Israel-Stewart defines the 
VAHydro scheme.^ vaHydro is expected to perform better than both IS theory 
and aHydro during all evolution stages. 


The vaHydro equations are obtained by generalizing the ansatz ( 1491 to include 


arbitrary (but small) corrections to the spheroidally deformed /rs( 2 ;,p): 


f{x,p) = fRs{x,p) -k 5f{x,p) = /i, 


^/Pf,Ef^^{x)p^ - fl{x) 
A{x) 


+ Sf{x,p). (156) 


The parameters A and fi in (156) are Landau-matched as before, i.e. by requir¬ 
ing {E)g = {E'^)g = Q. To fix the value of the deformation parameter ^ one de¬ 
mands that 6f does not contribute to the pressure anisotropy Pt—Pl] this requires 
{x^Xy+y^yu—2z^^Zi,){p^^p'''>)I = 0. Then, upon inserting (156) into (127), we obtain 
the vaHydro decomposition 




with 


subject to the constraints 

= {x^x^+yf,y^-2z^z^)n^'' = ff)) = 0 


(157) 


(158) 


(159) 


Clearly, the additional shear stress arising from 5f has only 4 degrees of freedom. 

The strategy in vaHydro is now to solve hydrodynamic equations for aHy- 
DRCP^ (which treat Pt—Pl nonperturbatively) with added source terms describing 
the residual viscous flows arising from 6f, together with IS-like “perturbative” equa¬ 
tions of motion for H, and The hydrodynamic equations are obtained by 
using the decomposition (157) in the conservation laws ([^. The evolution equations 
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for the dissipative flows H, and are derived by generalizing the DNMFpSl 
procedure to an expansion of the distribution function around the spheroidally de¬ 
formed /rs in (149), using the 14-moment approximation. These equations are 
lengthy; for massless systems undergoing (2-|-1)-dimensional expansion with longi¬ 
tudinal boost invariance they were derived by Bazow et Generalizations to 

massive systems and full (3-|-l)-dimensional expansion are in progress. We give 
their simplified form for (O-l-l)-d expansion in the next subsection. Especially at 
early times 5f is much smaller than 5f, since the largest part of 6f is already ac¬ 
counted for by the momentum deformation in (149). The inverse Reynolds number 
= -^TT^^'^TT^u/'Piso associated with the residual shear stress is therefore 
strongly reduced compared to the one associated with significantly improving 
the range of applicability of vaHydro relative to standard second-order viscous 
hydrodynamics. 


4.4. Testing different hydrodynamic approximations 


For (O-l-l)-d longitudinally boost-invariant expansion of a transversally homoge¬ 
neous system, the Boltzmann equation can be solved exactly in the relaxation time 
approximation (RTA), both for masslesP^Hini massive particles EIIlll More re¬ 
cently, an exact solution of this equation was also found for massless systems un¬ 
dergoing (l-l-l)-dimensional expansion]13111 with a boost invariant longitudinal and 
azimuthally symmetric transverse velocity profile discovered by Gubser (“Gubser 
flow” jMEU as the result of imposing a particular conformal symmetry (“Gubser 
symmetry”) on the flow. These exact solutions of the kinetic theory can be used 
to test various hydrodynamic approximation schemes, by imposing the symmetry 
of the exact solution also on the hydrodynamic solution, solving both with iden¬ 
tical initial conditions, and comparing the predictions of both approaches for the 
evolution of macroscopic observablesF^I^^I^^I^^I^^I^J 

We will here use the (O-l-l)-d case to test the VAHydro approach.l^ This il¬ 
lustrates the procedure and the kind of conclusions one can draw from such a 
comparison. Setting homogeneous initial conditions in r and space-time rapidity 
77s and zero transverse flow, reduces to a single non-vanishing component tt: 

= diag(0, —7f/2, —7f/2, if) at z = 0. The factorization nRs(C, A) = 7^o(C) ^iso(A) 
etc. are used to obtain equations of motion for A, if!^ 
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Piso(A). 
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Trei and the ratio of shear viscosity 77 to entropy density s, rj/s, are related by 
Trei = 5?7 /(sT) = 577/(7 ?.^/^(^)sA). The numerical solution of these equation^ can 
be compared with the exact solution of the Boltzmann equationf^ and also with the 
other hydrodynamic approximation schemes discussed above, plus a 3rd-order vis¬ 
cous hydrodynamic approximation!^ As an example, we show in Fig. [^the entropy 
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Fig. 1. (Color online) The particle production measure (T^n(Tj^))/(ron(ro)) — 1 as a function of 
Afnr\ls. The black points, red dashed line, blue dashed-dotted line, green dashed line, and purple 
dotted line correspond to the exact solution of the Boltzmann equation, vaHydro, aHydro, third- 
order viscous hydrodynamics,!^ and DNMR second-order viscous hydro dynamics p2l respectively. 
The initial conditions are ro = 600MeV, ^0 = 0, a-ud 'S'o = 0 at To = 0.25fm/c. The freeze-out 
temperature was taken to be Tj- = 150 MeV. 


production (measured by the increase in particle number rn(r)) between start and 
end of the dynamical evolution from an initial temperature of 600 MeV to a final one 
of 150 MeV. For this extreme (O-l-l)-d scenario, where the difference between longi¬ 
tudinal and transverse expansion rates is maximal, VAHydro is seen to reproduce 
the exact solution almost perfectly, dramatically outperforming all other hydrody¬ 
namic approximations. In particular, it should be noted that only the aHydro 
and vaHydro approximations are able to correctly reproduce entropy production 
(or rather the lack thereof) in both the extreme strong coupling (ideal fluid dy¬ 
namics, Treai = T^/s = 0 ) and the extreme weak coupling (free-streaming particles, no 
collisions, Treai, vls= —> oo) limits of the microscopic dynamics. vHydro schemes 
based on an expansion around a locally isotropic equilibrium distribution cannot 
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reproduce the constraint that entropy production should vanish as collisions cease; 
these schemes break down for large r]/s values. 

Similar comparisons have been done for the massive (0+l)-dimensional 
cas jUmUl] and for the (l+l)-dimensional Gubser flow. In all cases one finds the fol¬ 
lowing hierarchy of hydrodynamic approximations, when listed in order of improv¬ 
ing accuracy in their descriptions of the moments of the exactly known microscopic 
dynamics: first-order viscous hydrodynamics (Navier-Stokes theory), second-order 
Israel-Stewart theory, second-order DNMR theory, third-order viscous hydrodynam¬ 
ics a la Jaiswal, aHydro, and vaHydro. In view of the increasing sophistication of 
these approximation schemes, as discussed in the preceding subsections, this order¬ 
ing is not surprising, and some variant of vaHydro is likely to become the standard 
hydrodynamic modelling framework in the future. At the moment, however, only 
vHydro and aHydro have been implemented numerically for (2-|-l)-d and (3-|-l)-d 
expansion which do not rely on simplifying assumptions such as longitudinal boost- 
invariance and azimuthal symmetry. The fireballs created in heavy-ion collisions are 
not azimuthally symmetric, and experiments tell us that they feature characteris¬ 
tic anisotropic flow patterns that could never arise from an azimuthally symmetric 
initial condition. Longitudinal boost-invariance is not a good approximation ei¬ 
ther for particles emitted at large forward and backward rapidities, and it becomes 
worse when going to lower energies. Therefore, much effort is presently being ex¬ 
pended into developing (2-|-l)-d and (3-|-I)-d implementations of the aHydro and 
vaHydro schemes. 

5. Numerical Implementation of Hydrodynamics 
5.1. Need for t and rj 

The hydrodynamic simulations of the ultra-relativistic heavy ion collisions is best 
implemented in the hyperbolic coordinate systeni^ (also known as the Milne coor¬ 
dinate system) where instead of t (the laboratory time) and z (the beam direction), 
one uses the longitudinal proper time 


T = 

(161) 

and the space-time rapidity 


II 

1 

(162) 

Equivalently, t = rcoshry and z = rsinbry. 



One reason this is useful is that as shown in Fig. the evolution of the ultra- 
relativistic heavy ion collisions can occur only in the forward light cone bounded 
by the light cone axes, equivalently t = 0. More physically, suppose two identical 
systems were created at t = 0 and z = 0 by the initial collision of the two nuclei. 
Further suppose that one of the two has the collective velocity v in the z direction, 
but the other one is at rest. In this case, due to the time dilation, the same stage of 
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Fig. 2. A schematic diagram showing evolution of fireballs with differing Vz after the collision of 
two heavy ions at t = 0 and z = 0. 


the evolution will be reached when the lab time is at td for the system at rest and 
at tdj\/\ — for the moving system. In relativistic systems, these two Minkowski 
times can be very much different even though the two system are at the same stage 
in their respective evolution. However, in terms of the proper time both are at the 
same t = td since r is nothing but the local rest frame time. Hence, it is very 
natural that we should use the t — rj coordinate system when there is a strong 
longitudinal flow of matter. In the ultra-relativistic heavy ion collisions, this is the 
case due to the original beam momenta of the projectile and the target. 

For numerical implementation of hydrodynamics, one first needs to formulate 
the conservation laws in this coordinate system. For this, we need to know the 
transformation law between t — t] and t — z. We can start with the derivatives 


dr 


dt dz 

—dt + ^d^ 
or or 

cosh 7]dt + sinh rjdz 


(163) 


and 


dri = ^dt + ^dz 

OT] or] 

= T sinh ridt + r cosh ridz 


(164) 


which can be summarized as Lorentz transformations 


where da 


da = A^da and = A^da 


(9t, V_L, (1/t) 5,,) and 


/ cosh 77 0 0 — sinh r] \ 


0 10 0 

0 0 1 0 

\ — sinh ?7 0 0 cosh ry / 


(165) 


(166) 


The inverse transform is obtained by substituting t] with —rj. From now on, we 
will use the first letters of the roman alphabet (a, 6, • • •) to represent the components 
in the Milne space as above. 
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To figure out the transformation of the conservation law, we apply Eq.(165) to 


dtf + = cos\\T]drf - 


dr,f + dr,f - siuhridrf 


= - dr {t cosh T]j^ — TsinhTyj^) — dr, sinhryj° — - coshryj^) 


Defining 


and 


q'^ = cosh — sinh t] q^ 
(p = cosh qq^ — sinh q 


9“ = (9^,q-L,9’’), 
for any 4-vector the conservation law becomes 

0 = drirf) + Vj^-(rjj_) + dr,f = dairf) 
where only has the x, y components. One can also show 

Dr = U^d,r = U°-da 


(167) 

(168) 

(169) 

(170) 

(171) 


Note that our definition of the q component is different from the curvilinear defini¬ 
tion of the q component by a factor of r. We do this to keep the dimension of the q 
component the same as the other components. If one uses the curvilinear definition 

f = -r ( 172 ) 

r 

then the conservation law becomes 

0 = drirf) + Vj_-(Tj_L) + dr,iTf) = dairf) (173) 

where da = idr, V_l, 9,,). 

For the energy momentum conservation, one needs to apply the transformation 
law 3 times for the each index in dfT^'^ = 0. The algebra is a bit tedious but 
straightforwardlHUHHl The result is 


where we defined 


drT^'^ + -f d^T^'^ = -tr'" 


jab ^ 


and the index v = x,y. For the rq component, 

1 


1 


(174) 

(175) 

(176) 


drr^ + -f 

T r 

Both and conservation laws contain the geometrical source term. Trans¬ 
verse momentum conservation is simpler: 

= o 


(177) 
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where v and w are transverse coordinate indices. When testing one’s code for the 
energy-momentum conservation, the following form may be more convenient 


— / dxdydr] = 0 


(178) 


The half-transformed 7“'^ = satisfies the conservation law drT'^^ + drjT^^ + 

= 0 without the geometrical source terms exactly like Eq.(170|. 


For the shear and the bulk evolution equations in Eq.(145l or (147), transfor¬ 


mation from the Minkowski space to the r — ry coordinate space is straightforward, 
but it is a lot more involved algebraically. For more details see RefslHl^HH^ 


5.2. Numerical solution of conservation equations 

In this section, we will first discuss conservation laws in the Minkowski coordinates 
where the conservation laws are 




(179) 


In the Milne coordinate system, the energy-momentum conservation takes a slightly 
different form 




(180) 


where S° is the geometric source term in Eqs.(174| and (176). However, as will be 


shown shortly, the methods illustrated below can be easily adapted to this case. 

For vHydro, the energy-momentum tensor is given in a general reference frame 
by the decomposition 


rjiflU _ rpyU 




where 




(181) 


(182) 


is the ideal fluid part of the tensor. The net baryon current has the form 


Jb — Pbu^ + Vg. The equations that need to be solved are given in Eq.(179l, 


together with the relaxation equations for the dissipative flows, for example the 


Israel-Stewart equations (147). The first step in solving the hydrodynamic equa¬ 


tions is their initialization. Let us assume that some microscopic pre-equilibrium 
dynamical theory provides us with a baryon current Jb{^) 3-nd energy-momentum 
tensor for points x^ on some Cauchy surface on which we want to initialize 

the hydrodynamic evolution stage. The following projection steps, to be taken at 
each point x on that surface, yield the required initial value fields^ 

First we define the local fluid rest frame by solving the eigenvalue equation 
= eu^ for its normalized timelike eigenvector . The associated eigenvalue 
gives us the LRF energy density e. The LRF baryon density is obtained from Jg 

‘The following paragraph refers to the vHydro decomposition | |l0|l8l[ |. A slightly modified pro¬ 
jection method applies for the aHydro and vaHydro decompositions 11501, 11511 and (157 
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by projecting onto u^\ ps = u^J^. The initial heat flow Vg is the component of Jg 
perpendicular to u^: Now that we know the LRF energy and baryon 

densities we can compute the LRF pressure P from the equation of state of the 
fluid, P{e, pb)- Next, the bulk viscous pressure is obtained from If = — ^ — 

P. Finally, the shear stress is obtained as — eu^u'^ + (P+n)A^'^ or, 

equivalently, as = A^^T“^. 

For a simple illustration of numerical methods that can solve the hydrodynamic 
equations, let us first consider a single conservation law in 1-D. There is no difficult 
conceptual obstacle in extending this case to the multi-component, multi-dimension 
cases such as the Israel-Stewart equations. The conservation equation is 

dtu = -dooj (183) 


We need to supplement this equation with a relationship between the density u and 
the current j. Simplest example is j = vu with a constant speed of propagation 
V. But in general j is a non-linear function(al) of u. For instance, the ideal part, 
Eq.(182), is certainly not in this simple form due to the normalization condition 
Uq = Vl — and also to the presence of the pressure term. In the dissipative 
cases, the relaxation equation 


{dt -f l/TR)j = - (D/TR)dxU (184) 

determines the relationship between j and u. In such cases, the numerical methods 
discussed in this section needs to be applied in two steps. In the first step, the 
conservation laws are used to advance the time component of the currents using 
the methods that will be discussed here. In the second step, the spatial part of 
the currents needs to be reconstructed from the time components. The relaxation 
equations also need to be solved separately, although the techniques discussed here 
can be easily adapted to handle the relaxation equation as well. 

For the simple j = vu case with a constant v, the equation becomes 

dtu = —vdxU (185) 


This is an advection equation and has a simple solution 


u{t,x) = f{x — vt) (186) 

That is, at any given time, the solution is just the translation of the initial profile 
by vt. Analytically, this is trivial. However, it is remarkable how difficult it can be 
to maintain the initial proHle in numerical solutions. We will often use this as the 
simplest test case for our algorithms. 

To solve the conservation equation numerically, one first needs to discretize time 
and space. We define 


u2 = u{tn,Xi) (187) 

where = to + nh and Xt = xq + ia for any function u{t, x) of time and space. 
Here h is the time step size and a is the spatial cell size. Physically, it is important 
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to have the following properties in the discretization method. First, the method 
should conserve total u explicitly, that is, J2i "^7 — J2i modulo the boundary 
terms. For this, one requires that the discretized form of the divergence to take the 
form 


idj)l 




(188) 


where ji{u"') is the discretized representation of the current j at Xi ant One 
can easily see that in the sum J2i=oi^xj)i only the boundary terms would survive. 
The details of the boundary terms depend on the method. However, as long as the 
boundaries are far away from the physical region, the boundary terms should be 
vanishingly small. If the boundary of the space is not too far away from the physical 
region, then some suitable discrete boundary conditions should be imposed. 

The second requirement is simple, yet quite demanding: If u® > 0 for all i, then 
we would like this property to be maintained for any future time. For instance, if 
u represents the energy density, then it should never become negative. 

To illustrate some of these issues, consider again the advection equation dtu = 
—vdxU with u > 0. Let the initial condition be a rectangle: = Uc for b < i < f 

and = 0 otherwise. This is a prototype of many situations where two smoothly 
varying regions are joined by a stiff gradient. The simplest discretization method 
of dtU = —vdxU is the forward-time centered-space (FTCS) method 


n+1 


= Ui - 


vh 

2a 




(189) 


which is correct up to 0{h?) and 0{a^) errors. Since the second term in the right 
hand side is supposed to be a correction, we require \vh/a\ < 1. This is certainly 
in the form of Eq.(188) and hence conserves the total u. 

According to the analytic solution Eq.( 186), the space behind the back edge (xt) 
of the rectangle should always have u" = 0. However, according to Eq.(189) the 
cell right behind the back edge becomes non-zero and negative after the first time 
step 


Ub-l = 


vh 

-1 

2o 


(190) 


At the same time, the front edge starts to get distorted by the same amount 


uj = Uc- ul^ > Uc 


At the next time step, u'^_i becomes even more negative 


2 vh 

Ub-l = 


u/i Y 
^ ) 


Uc 


while Uf deviates even more from Uc 

2 2 
Uf = Uc - Ub_i 


(191) 


(192) 


( 193 ) 
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Furthermore at t 2 , the next cell in the empty region becomes non-zero = 

Uc- Clearly, both the profile preservation and the positivity of the solution 
are grossly violated by this method even though the total u is conserved. In addition, 
one can easily see that the growth at Xb-i will continue, indicating that this method 
is unconditionally unstable. 

In order to cure the negativity problem, one may note that the trouble above 
mainly comes from the centered nature of the O(a^) numerical derivative in 
Eq.(189). Instead, one may try the first order approximatior0 


n+1 


vh , 


a 


(194) 


In this “up-wind scheme”, ul_^ trivially vanishes at ti. In fact, all Ub-i for z > 1 
will remain zero for all times. Similarly, one can show that the numerical solution 
is positive and bounded everywhere as long as \vh/a\ <1. If ?; < 0, mirror-image 
conclusions can be reached if one uses {dxu)i « (ui+i — Ui)/a. 

Positivity is thus maintained in this up-wind scheme. However, as the system 
evolves in time, the shape of the solution gets more and more distorted. This is 
because the first order difference 


« - = dxu{t„,Xi) - -dlu{tn,Xi) -b O(a^) 


(195) 


is too crude an approximation of the first order derivative. The second derivative 
term in Eq. (1951 in fact introduces too much artihcial (numerical) damping to 
preserve the shape for long. In effect, the profile at t is the convolution of the initial 
profile and the Gaussian Green function of the diffusion equation (Eq.(92)) with 
the diffusion constant given hy D = ^. As one can see in Eq.(92), the width of the 
Gaussian grows linearly with time. Hence, the initial profile will be smeared out 
quickly. 

Better discretization methods must keep the positivity preserving nature of the 
up-wind method and at the same time must have a better approximation for the 
spatial derivative than the simple first order difference for shape preservation. To 
devise better discretization methods more systematically, consider hrst dividing the 
space into N + 1 cells of size a labelled by integers 0 through N. The z-th cell starts 
at Xi-i /2 = Xi — al2 and ends at Xi+ 1/2 = Xi + a/2. (See Fig. Averaging over 
one spatial cell and integrating over one time step, the conservation law dtu = —dxj 
becomes 


-n+1 


1 


dt 


■,^1+1/2) ~ Xi_i/2)^ 


(196) 


where u/ = Ui{tn) is the cell average. This is an exact expression which can be 
used as the basis for further approximation. This type of methods are known as 
the finite volume methods. 

j This is in fact trivially exact when a = vh since in that case That is, the whole 

profile moves by one spatial cell at each time step. However, this wouldn’t work for more general 
currents. 
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j(i-l/ 2 ) j(i+l/ 2 ) 



Fig. 3. Spatial grid used in deriving finite volume methods. The solid line is the lowest order 
histogram function approximation of u{t,x) in the x^-centered grid. The dot-dash line is the 
projection of the x^-centered grid onto the grid. 


Approximating the time integral using the midpoint rule, we get 

= w" — — ^j(in+l/2: a^i+ 1 / 2 ) — jVnH-1/2, +0{h^) (197) 


where tn+i /2 = tn + h/2. There are few things that should be mentioned here. 
First, the basic quantities to calculate are the cell-averaged values it". Second, we 
need to approximate the function it(t, x) itself from u" because the right hand side 
contains u{t,x) evaluated at t„+i /2 and Xi±i/ 2 - An approximate form of u(tn,x) is 
also needed for obtaining for the next time step (see below). Therefore, how 
we approximate u{t, x) and evaluate determine different numerical schemes 

and the accuracy of the given scheme. 

In many schemes, the values at spatial half points are not unique because the 
interpolating function approximating u{t^ x) is usually only piece-wise continuous. 
For instance, see Fig.j^and Fig.|^ One way to deal with the ambiguity in evaluating 
j(t„+i/ 2 , a^i±i/ 2 ) is to just avoid evaluating it at the boundaries. The staggered 
nature of the space and time indices in Eq.(197) suggests an obvious way to do so. 
Suppose that the initial data is given for u" where the i-th. grid is centered at Xi and 
has the spatial interval 1 / 2 , 2 ^ 1 + 172 ] • For the next time step, instead of updating 
the values of Ui within the intervals [xi_if 2 ,Xi+i/ 2 ]i we update the values of Mi+i /2 
within the shifted intervals [xi,Xi^i]. Using the lowest order approximation for 
u{t,x) (the histogram functions in Fig. [^, this yields 

‘ h 




u 


^i-\-l 


1+1/2 




(198) 


where (m” -I- Mf_^]^)/2 is the average value of u{tn, x) in the shifted interval [xi, Xi+i] 
using the histogram function in Fig. 3 For the values of at tn+i/ 2 , the 
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forward Euler method 




(199) 


is enough since the overall error in the midpoint rule is 0(h^). In this way, we 
avoid evaluating the current at the cell boundaries. For the next time step, we shift 
back to the original grid. This staggered method is a slightly generalized form of 
the Lax-Friedrichs scheme which is second order accurate in space and third order 
accurate in time. The original Lax-Friedrichs scheme is only second order in time 
because the currents are evaluated at instead of at t„+i/ 2 - 

Eq.(199) involves evaluating numerical derivative of the current j at Xi. This 
cannot be exactly determined when all one has is data on the average values ft". For 
instance, {jx)2 could be the first order approximations (j7+i ~ jf)/a or (jT-j”_i)/a, 
or the second order approximation — j^_i)/2a, or any other approximate form. 
In normal situations, choosing a higher order formula should be better than the first 
order ones. But this is not always the case when the gradient is stiff. 

We have already shown that using the central difference ~ in the 

forward Euler method (the forward-time-centered-space method in Eq.(189l) can 
be disastrous when the gradient is stiff, although it can be safely (and preferably) 
used in smooth regions. When gradient is stiff, one should instead use the up-wind 
method Eq.(I94) to maintain the positivity. Therefore, an intelligent scheme would 
choose the derivative according to some approximate measure of the true gradient. 
One way to do this is to choose the gradient according to the following scheme 


(5.w)r = 


0if<<12r±i 


> u: 


i±l 


else sign(M”_^;^ — m”) min(0- 




( 200 ) 


The first line indicates that the function has either a maximum or a minimum within 
the interval. Therefore the slope at Xi can be best approximated by 0. The second 
line applies when the function is changing monotonically near Xi. The parameter 
1 < 0 < 2 is there to be slightly more general. This choice of the derivative is called 
the “generalized minmod flux limiter”. 

The Lax-Friedrichs scheme represented by Eq.(I98) is only 0{a?) accurate be¬ 
cause we have used the histogram function as an approximation for u{tn,x). Need¬ 
less to say, this is the lowest order approximation. As a result, the Lax-Friedrichs 
scheme contains too much numerical diffusion to be practically useful. Both of these 
facts can be easily seen from the Taylor expansion 


H+l 


A-HI/2 g 


-rrd: 


;<+i/2 + 0{a^) 


( 201 ) 


where the second derivative term is the 0{a^) error term that also causes strong 
diffusion. In time, this diffusion distorts the solution too much just as in the first 
order up-wind method. 

To obtain a better approximation, one needs to evaluate u{t, x) more accurately 
using the average values from the nearest neighbor cells. If one uses the linear 
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x(i-l/2)+c(i-l/2)dt 



Fig. 4. A schematic view of the cell division used in the Kurganov-Tadmor scheme. 


interpolant (see Fig. 13 


u^{x) = It” + {dxu)^{x - Xi) 


( 202 ) 


for each interval 2 ;i_|_i/ 2 ], then the scheme should become at least 0(a^). 

Using this to improve the estimate of w(Yi /2 ^^ids a correction term to the lowest 
order result Eq.(198) 


f.n+l 

“i+l/2 


^i+l 


h 


i(tn+l/2yXi+l) ~ i/(^n+l/2) 


- g ((a.u)r+i - {d,u)-) 


(203) 


This staggered approach is the second order NT (Nessyahu-Tadmor) scheme and 
works reasonably well.l^ The last term in the right hand side of Eq.(203l cancels 


the second derivative term in Eq.(201| so that together they represent it'Yi /2 
the O(a^) error. In other words, the last term in Eq.(203l is the anti-diffusion term 


that is correcting the large second order diffusion introduced by the symmetric 
combination (u" -|- -u'Yi)/2- Therefore, this scheme is O(a^) in the smooth region. 

Why don’t we then just replace these terms with If one does that, it 

just becomes the forward-time-centered-space (ETCS) scheme in Eq.(189|. The 


difference, the 0{a'^d^u{x)) term, provides just enough numerical diffusion so that 
spurious oscillations do not propagate from the stiff gradient. 
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The second order NT algorithm represented by Eq.(203) is related to another 


Eq.(203) but replace 


often-used finite element method called the SHASTA algorithm. We start with 

r+i/2 + o (^*+1 ~ -I- Ui) + 0{a^) (204) 


f.b 

“i+1 


SO that 


^”+1% - 


i+l/2 ^(j{tn+l/2T j{tn+l/2t 

- 2<+i/ 2 + O - ^ ((9.«)r+i - 


(205) 


At this point, this is no longer a finite volume method even though we will keep 
using the notation uf to represent the value of u at and Xi. It is also not a 
staggered method any more as the right hand side contains 'S(Yi/ 2 - Re-labelling the 
spaces indices i + 1/2 —i’ i and i —)■ i — 1 so that the new grid size is a' = a/2, we 
get 

^/(^n-l-1/2 ; ^j-l-1) J (^n-|-l/2) 1 

+ g(zZ(Vi - 2u/ + uU) - ^ ((9,u)r+i - {d,u)U) (206) 

with the appropriate scaling of the derivatives in the last term and after renaming 
a' —)■ a. With u(t„_|_i/ 2 ,Xi±i) given by Eq.(199), Eq.(206) represents the basic 
SHASTA algorithm. 


In practice, Eq. (2061 is broken up into two stages to ensure positivity. Specifying 
j = vu, the first transport step i^ 


= u/ 


1 


[u: 


- ^(<+1 - Y(9.«)r+i - uu + Y(5x«)r_i) 


i+l + U^_i) 

vh, 

J 


vh , 

T 


(207) 


again using Eq.(199|. If u/ > 0 for all I, then as long as \vh/a\ < 1/4, 


,,n+l 


IS 


positive. 

The second stage is the anti-diffusion step 






(208) 


where represents numerical estimate of the second derivative at Xi that 

preserves the positivity. The numerical approximation suggested by the last term 


in Eq. (2061 turned out not to preserve the positivity. The original formulation of 

(209) 


the SHASTA algorithm by Boris and Book uses a conservative form 

1 


idiw)r^ = - - iw,)r^) 


^ In the original SHASTA algorithm, is used in place of 

{^dx'a)'2j^^ — /2 in this stage. In practice, as long as \vh/a\ is small, this difference 

does not matter much. But it is crucial that a flux limiter is used in the second anti-diffusion step. 
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where 


Mti = 


0 if and do not all have the same sign, 

else sign(Ar+i)min(8|Ari|, |Af^^^8|A^^^|) 


( 210 ) 


with 


A"+i = «+i-<+i)/a 


( 211 ) 


This is similar to the minmod flux limiter and maintains the positivity. The 
SHASTA algorithm is used in many hydrodynamics simulation programs for ultra- 
relativistic heavy ion collision^^ including pioneering works in RefslHUH®! and also 
later works in RefsE3HlHl2l] 

The second order NT scheme and the SHASTA scheme in practice work fairly 
well 11221 However, in these schemes the ft. —>■ 0 limit cannot be taken since the 
numerical viscosity behaves like ^ 1/ft. It would be convenient to be able to take this 
limit because one can then formulate the discretized problem as a system of coupled 
ordinary differential equations in time. Many techniques for the ordinary differential 
equations such as the Runge-Kutta methods then become available to control the 
accuracy of the time evolution. So far different time integration techniques were 


not available other than the midpoint rule in Eq.(197). 


One way to achieve this is to subdivide the cells as shown in Fig.j^with the piece- 
wise linear approximation for u{t, x). The size of the cell containing the discontinuity 
at the half integer point Xij^i /2 is controlled by the local propagation speed Ci+i/ 2 - 
That is, the subcell surrounding Xi +112 is between the left boundary = 

Xi+ 1/2 — Ci+i/ 2 ft and the right boundary = Xi+ 1/2 -f Ci+i/ 2 ft- Then the cells 

containing the boundaries and the cells not containing the boundaries (between 
^i-i /2 ^i+ 1 / 2 ) independently evolved. For the subcell containing Xi+ 1/2 




n+l 

"i-l-1/2 


i-l-1/2 + '*^i+l/2 


2c. 


■ 1 + 1/2 


‘‘i+lj'l) j(^i+l/2) ) + ^(^) (212) 


which is basically Fq. (1981. Here the superscripts r and I means the value of the 


approximate M(t„, x) at the boundary points a;/_|_ 2/2 ^i+ 1 / 2 ’ respectively. Within 

the smooth region between ^i+ 1 / 2 ’ using Eq.(197) 


,,n+l 


“i+l/2 + ^i-1/2 


(j(1iUi/2)-J(«[-i/ 2)) +0(ft") (213) 


where the first term in the right hand side is obtained by applying the trapezoid 
rule. The divided cells are then projected onto the original grid using the size of 
the cells as the weight to get 






_ (c»-l/2 + Ci+l/2)ft \ „+i 
a ; 


( 214 ) 
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When the h —>■ 0 limit is taken, this procedure yields 

dUi _ -ffi+l/2 ~ 


dt 


where 




i+l/2 


^■(w+i/2)+j(u,+ i/2) 


- C, 


t+1/2- 


U'^, 1 /o — M - , 1 /o 
2 + 1/2 2 + 1/2 


(215) 


(216) 


where now values of the piece-wise linear u{t,x) when 

approaching Xi^i /2 from the right and from the left, respectively. They are given 

by 


,-7 + 


— '^i+1 (c^/2) (<9a:ri)2-j-l 

^2-t-l/2 — 


(217) 

(218) 


again using the minmod flux limiter for the derivatives. This is known as the second 
order Kurganov-Tadmor (KT) schem^iSSl and it is implemented in the 3-l-lD event- 
by-event viscous hydrodynamics simulation program Musicl21221 The numerical 
viscosity in the smooth regions is known to be O(a^d^u). 

The structure of the KT algorithm, Eqs.( |215] ) and ( |216 ), is the same for the 
lowest order, second order and the third order algorithms. All one needs to do to 
improve accuracy is to get a better estimate of Actually, it is instructive to 

consider the lowest order result which uses the histogram function as the approxi¬ 
mation of u{x). One then has “ ^i+i m~|-i / 2 = Ui. It is easy to see in 

this case that if j = vu and hence Ci_|_i /2 = |'a|, Eq.(2151 automatically becomes the 
up-wind method. 

We now have a set of ODE’s. What is a good choice of the ODE solver? One of 
the physical requirement is again the positivity. For instance, suppose u represents 
particle density. One knows that u can never be negative. One of the schemes 
that preserves the positivity is the Heun’s method which is one of the second order 
Runge-Kutta schemes. The equation du/dt = / is numerically solved by following 
these steps 


u* =u'l + hf{tn,u'^) 


u** =U* + hf{tr,+uU*) 

-n+l 


= 


(219) 

( 220 ) 

( 221 ) 


Positivity is maintained at each stage with a suitable choice of the flux limiter such 
as the minmod flux limiter. 

For 3-D (and similarly for 2-D), a simple extension 


d _ _ ^i+l/2,j,k ^i-l/2,j,k 

~T7^iik — 

at a-r 


try _ zjx 

-^2,j + l/2,fe ■‘^i,j-ll2,k 


TTZ _ TLTZ 

^i,j,k-\-l/2 '^2,J,fc-l/2 


( 222 ) 
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works well. In curvilinear coordinate systems such as the Milne coordinates, one 
may not have the conservation law in the form = 0. Instead, it may take the 

form 


dar = s 


(223) 


where S is the geometrical source term that does not involve derivatives. The KT 
algorithm can be easily adapted to this case by simply adding the source term on 
the right hand side. Namely, 


dui 

dt 


Ht+l/2 - 

-h S{Ui) 


(224) 


In implementing the KT scheme, one needs the maximum speed of propagation 
at Xi±i/ 2 - If there is only one variable u and one current j, then the speed of 
propagation in the i-th direction is 


du 


(225) 


If there are more than one current, then we first define the Jacobian matrix 


jab ^ ^ 

* du. 


(226) 


where i = 1, 2, 3 is the space index and a = 1, • • • , M labels the conserved quantities. 
Therefore we have 3 M x M matrices. The maximum propagation speed in the i-th 
direction is 


Ci = max(|Ai|,-• • , IAa^I) (227) 

where A’s are the eigenvalues of the Jacobian J“^. When discretizing in time, the 
original authors of the KT algorithm recommended the time step h to be small 
enough so that |max(ci)/i/a| < 1/8. 

For the explicit form of the Ci for the 3+ld ideal hydrodynamics including the 
net baryon currently used in MusiC, see Ref.^^ For implementation of the event- 
by-event 3-|-ld viscous hydrodynamics in Music, see RefsEEH Additional valuable 
information on the algorithms used in solving the vHydro equations can be found 
in Ref.ESl For some comparisons of various schemes discussed in the section, see 
Ref.ES^ 

There are many other numerical schemes that are currently in use but we un¬ 
fortunately did not have space to discuss. These include, but not limited to, a 
PPM (Piece-wise Parabolic Method) scheme,^^^ a Lagrangian scheme where the 
grid points follow the movement of fluid cells,^^^ Riemann solversand a SPH 
(smoothed particle hydrodynamics) method.l^^ The simple FTCS scheme has also 
been employed for the viscous hydrodynamic£2lIIi3 for smooth initial conditions. 
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Fig. 5. Plot of T(r, r)/To for the Gubser solution. Here qvQ = (3 + '\/5)/2. The left panel shows 
T(r, r)/To at fixed r’s and the right panel shows T(t, r) at fixed r’s. 



Fig. 6. Plot of the hypersurface with T(t, r)/To = 0.17 for the Gubser solution. 


5.3. Freeze-out Hypersurface and Cooper-Frye Formula for Particle 
Production 

When a hydrodynamic system is expanding, there comes a time when the system is 


too dilute to be treated with hydrodynamics (c.f. section 4.1). From this point on, 


the system is basically a collection of non-interacting particles. In realistic systems, 
this time is not the same for all fluid cells. As the system starts to expand at tq, 
there are cells at the periphery of the system that are dilute enough to “freezes 
out” in a very short time. As time elapses, expansion of the system can cause 
hot and dense matter to flow into the location of those frozen cells. These will 
eventually freeze-out, too. Therefore, the freeze-out surface volume cannot be a 
simple 3-dimensional volume. It is a complicated 3-d volume in the 4-d space-time. 

To illustrate this point, consider the following Gubser solution of the boost- 
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invariant and azimuthally symmetric ideal- hydrodynamics!^^ 

2To 


r(r,r) = 


(2qT [1 -I- + r^) -|- — r^)^] ) 


1/3 


(228) 


where = x^+y^ and r = — z^. It is a simple matter of taking the r derivative 

of T{T,r) to see that for r > tq where rg = there are two values of r where 

drT vanishes. Since drT is negative for small r, this means that T{T,r) at fixed 
r > rg will have a minimum and then a maximum. This is illustrated in the left 
panel in Fig. The solid line is for qr = 1 which is near the center of the system. 
The temperature at that position decreases monotonically. At r = tq, there is an 
inflection point but the behavior is till monotonic. In the qr = 5 case, one can 
clearly see that the temperature decreases at first but a some point it starts to 
rise again as the pressure pushes hot matter from the central region towards the 
periphery of the system. Assuming that the freeze-out temperature is between the 
minimum and the maximum of T/Tg, the position qr = 5 will contribute to the 
final particle spectrum (freeze-out) three times. 

The plot of isothermal curve with T{t, r) jTg — 0.17 is shown in Fig. One may 
think of this as representing the “freeze-out” hypersurface in the Gubser solution. 
The long tail that can be seen above gr > 5 is unrealistic. In more realistic simula¬ 
tions, the simulation starts at times above the long tail. Nevertheless, the Gubser 
solution contains much of the features that the more realistic numerical solutions 

exhibitEEnitnii 

At the freeze-out hypersurface, the fluid is to be converted to particles according 
to the local equilibrium condition. In the kinetic theory, the particle number current 
for the i-th particle is given by 

if (a;) = ft y h[x,p)p>^ (229) 


where fi{x,p) is the phase space density and gi is the degeneracy. The number of 
Tth particle in a 3-d hypersurface is then given by 



(230) 


where da^ is the hypersurface element. Hence, the momentum spectrum is 

^ J ( 231 ) 

This is the celebrated Gooper-Frye formula.!^^^ 

Hydrodynamics deals with the energy density and its flow. Experiments measure 
the momentum distribution of identified particles. Gooper-Frye formula provides 
the link between them. In ideal hydrodynamics, local equilibrium is strictly main¬ 
tained. Hence once we find the freeze-out hypersurface (equivalently the energy 

* Analytic solutions of the ideal hydrodynamics do exist for special cases, but not for general cases. 
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density via the equation of state), all one has to do is to integrate Eq.(231| over 
the freeze-out hypersurface with fi{x,p) = -|- at) where Tpo is 

the freeze-out temperature, Oi = — 1 for bosonic statistics, Oi = 1 for fermionic 
statistics, and = 0 for classical (Boltzmann) statistics. 

The freeze-out surface we need to integrate over has the shape shown in Fig. 
which requires a closer inspection. Suppose all parts of the system reaches the 
freeze-out temperature at the same Minkowski time tpo and only once. In that case, 
the freeze-out hypersurface is just the Minkowski spatial volume. The integration 
element is just then dao = dxdydz. All other components vanish. This is simple, 
but not a realistic scenario in relativistic heavy ion collisions as explained in the 
previous section. 

In the Milne coordinate system, the volume elements in each direction are given 

by 


daa = {rdxdydri, —Tdrdydr], —Tdrdxdrj, —rdrdxdy) 


(232) 


When the freeze-out hypersurface is specified by the freeze-out proper time 


Tf(x, y, rj), the surface element on this surface is obtained by replacing dr in Eq.(232) 
with [dxf /dxi)dxi 


dal = ( 1 


dxf dxf dxf 
dx' dy ' dr] 


Tfdxdydr] 


(233) 


If dal is time-like, then it is guaranteed that p°'dal > 0. Therefore, Eq.(231) has 


a well defined interpretation that the particles are flowing out of the hypersur¬ 
face. If dal is space like, then depending on the direction of p“, p°‘dal could be 
either positive or negative. Equivalently, particles can be flowing into or out of the 
hypersurface. 

Before we discuss the physical situation when this can happen, first we define 
the dynamic rapidity y = [pz/E) so that 


E = rriT cosh y 
p^ = rriT sinh y 

with mp = The Milne space momentum is 

= (mrCOsh(?/- yy), pp, tot sinh(?/- ?7)/r) 


(234) 

(235) 

(236) 


which is equivalent to using from (c.f. Eq.( 166)) but with an extra factor of 1 /r 


for the T] component to conform with the geometrical definition of the hypersurface 


element Eq.(233). 


To illustrate what happens when dal is space-like, let’s consider the Gubser 
solution again. In this case, r/ (r) is a function of the transverse radius r = 
only. Since the system is boost-invariant, we can also set 77 = 0 without loss of 
generality. The inner product of the momentum and the volume element is 

p'^dal = {rriT coshy — {pT-h)drTf) Tfdxdydr] 


(237) 
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where f is the unit vector in the transverse direction. This can become negative if 
dTf/dr > {rriT/ pt) coshy > 1. From Figs. andone can see that large positive 
gradient drf /dr exists in the region where the temperature reaches Tpo the second 
time. At this time, the temperature is going up from the minimum in the left panel 
of Fig. Therefore, this part of the Tf is not about the fluid freezing-out. It is 
rather colder matter being heated up to become the a part of the fluid again. 

Before one could use the Cooper-Frye formula, one needs to know the freeze-out 
hypersurface. This is not a trivial problem because hydrodynamic simulations only 
provide the freeze-out space-time points. The freeze-out hypersurface needs to be 
reconstructed from these points. In 2-|-ld, this is not so complicated. But in 3-|-ld, 
it can become very involved. Discussion of this important topic however is beyond 
the scope of this review. An interested reader should consult Ref.l^i^ 


How do we use the Cooper-Frye formula Eq.(231|? If hydrodynamics is not 


coupled to the hadronic cascade after-burner, then usually the hypersurface inte¬ 


gral in Eq.(231) is performed as it is after the hypersurface is reconstructed. The 


rationale behind it is that when the cell crosses the freeze-out boundary 3 times, 
the contributions from the first two times should largely cancel each other. This 
is physically sound because when the fluid cell crosses the freeze-out surface the 
second time, the particles that are being heated up again are the remnants of the 
first crossing. When coupling to the hadronic after-burner, an additional condition 
such as p°‘da[ > 0 is usually employed .IHSHllSI 

The presence of non-zero shear tensor and the bulk pressure H signals non¬ 
equilibrium. In this case, the Cooper-Erye formula needs to be modified to take 
into account the non-equilibrium phase space density. Let 

f{x,p) = feq{x,p)-\-6f{x,p) (238) 

Then the consistency between the hydrodynamic and the kinetic theory = 


/ ( 2 ^) 3 po Py/ requires 


(Pp 


(27r)3p< 


-pinp'') Sf{x,p) 


and 


n = - 


m 


(Pp 


(27r)3p( 


df{x,p) 


(239) 


(240) 


with 


= K^pp'^p^ (241) 

These conditions can be satisfied by 

dfix,p) = (^A{Ep)n{x,p) -\- B{Ep)p''>^p‘'Kp^{x,p) H-) feq{x,p){l a^feq{x,p)) 

(242) 
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where again Ep = p^Up, is the energy in the local rest frame. Here A{Ep) and B{Ep) 
must be consistent with Eqs.(239l and (240), but otherwise arbitrary at this point. 
In the presence of the net-baryon current Eq.(242) also includes C{Ep)p ^'^'>with 
the constraint on C{Ep) given by Eq.(138). Exact forms of A{Ep) and B{Ep) (and 
C{Ep)) depend on the form of the scattering cross-sections in the underlying micro¬ 
scopic system. In Ref.it is argued that for most theories, the Ep dependence of 
A and B should be between 1 and 1/Ep. See also ptefsilMiHU In most simulations, 
the constant ansatz is used. 


6. Summary 

In this review, we have tried to deliver a more general and pedagogic view of the 
relativistic hydrodynamics currently used in the study of ultra-relativistic heavy ion 
collisions. One message we tried to convey to the reader was that hydrodynamics is 
a very general framework and yet it can describe a vast set of complex phenomena. 
Especially in QGP studies, hydrodynamics is an indispensable tool that connects 
the QGP properties to the actual observables. 

Another message we tried to convey was that the theory of hydrodynamics is 
a fascinating subject by itself. As discussed in section and section there are 
still many unsolved problems such as finding Kubo formulas for the second order 
transport coefficients and finding the right anisotropic equation of state. In view of 
the apparent collectivity in the high-multiplicity proton-proton and proton-nucleus 
collisions at the LHG, the theory of collective motion in small systems is also in 
urgent need of development. In these systems, thermal fluctuations during the 
hydrodynamic evolution may not be completely ignored ]l22i\12§i hope that this 
review has provided interested readers enough starting points to pursue these and 
other interesting topics on their own. 

In any short review, omission of some important subjects inevitably occurs due 
to the lack of space. One notable omission in this review is the discussion of the 
initial state models. As briefly discussed in section the initial condition of the 
hydrodynamic evolution must be given outside of the theory of hydrodynamics. 
For a meaningful comparison with the experimental data, having the right initial 
condition including the right energy-momentum fluctuation spectra is crucial. Un¬ 
fortunately, it is outside the scope of this review and must be left as the subject of 
a future review. 
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